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This dissertation examines the effects of ocean mesoscale variability on acoustic 
arrival time patterns for two separate ocean environments. First, for an open ocean 
environment away from strong boundary currents, the effects of randomly phased linear 
baroclinic Rossby waves on acoustic travel time are shown to produce a variable overall 
spreading in the arrival pattern, primarily producing a delay in the later, axial arrivals. 
Second, using the state-of-the-art Semtner-Chervin eddy resolving global ocean circu- 
lation model coupled with the University of Miami Parabolic Equation (UMPE) acoustic 
propagation model, the effects of a fluctuating frontal region created by the California 
Current on the temporal, spatial and seasonal variability in the individual modal arrivals of 
the first thirty modes over a one-model-year time span is assessed. The mesoscale bias 
variability is also examined by comparing the various peak arrival times for the range- 
averaged environment to that of the range-dependent environment. To support this work, 
approximate “wide angle PE mode functions” were newly developed which form a 
different basis set for modal expansion from that obtained using standard normal mode 
theory. These new mode functions provide the proper basis set for modal expansion of the 
field computed by wide angle PE models. 
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I. INTRODUCTION 



Knowledge of ocean acoustics and associated environmental parameters is 
essential to developing and strengthening the Navy’s ocean environmental/acoustic 
nowcasting and forecasting capabilities. As naval systems, platforms and sensors become 
increasingly more sophisticated, knowledge of environmental factors and their impact on 
acoustic propagation will play an even more critical role in undersea warfare (USW) 
system performance. To understand, predict, and assess the performance of various USW 
acoustic systems, it becomes essential to understand the temporal and spatial variability of 
the sound speed structure, the impact of oceanographic features such as fronts and 
mesoscale eddies on acoustic propagation, and acoustic pulse arrival time variations due 
to mesoscale features, internal waves, bottom properties and ambient noise. This disser- 
tation focuses specifically on the effects of ocean mesoscale variability on acoustic arrival 
time patterns. 

Oceanic mesoscale phenomena are the analogues of weather systems. Mesoscale 
features extend on the order of 100-250 km and persist on the order of 50-150 days. More 
than ninety percent of the kinetic energy of the entire ocean may be characterized by 
mesoscale variability rather than by large scale currents (Spindel and Worcester, 1990). 
These mesoscale features can significantly influence acoustic propagation and thus impact 
tactical operations. The impact of mesoscale features on the acoustic field at the receiver 
can be significant, affecting the temporal sequence, amplitude and angular distribution of 
the acoustic energy. 

The presence of mesoscale phenomena was clearly established by oceanographic 
studies such as the MODE expedition in 1973 (Mode Group, 1978). A method for 
monitoring this mesoscale variability in the ocean interior, ocean acoustic tomography, 
was introduced by Munk and Wunsch (1979). In this method, the sound-speed fluctuation 
field within a volume is measured by the travel times of variously oriented acoustic trans- 
missions through the volume. The acoustic travel time is inversely proportional to 
average sound speed along the transmission path. Shallow mesoscale perturbations will 



l 



alter only the average sound speed for steep ray arrivals, whereas mesoscale perturbations 
which extend to sound speed axis depths will alter the average sound speed for both steep 
and axial rays. Thus, by measuring travel-time perturbations of acoustic signals along 
various paths, the underlying sound speed structure may be inferred by inversion. The 
first at-sea experiment of this technique was performed in 1981 (Ocean Tomography 
Group, 1982). Subsequent efforts to refine tomographic techniques led to a number of 
investigations into travel time perturbations. A review of these studies is provided as 
background in Section A to help put into context the contributions presented by this 
dissertation. 

Understanding inherent ocean variability is also important when installing 
"permanent" acoustic transmission systems, such as those considered in support of the 
Acoustic Thermometry of the Ocean Climate (ATOC) (Munk and Forbes, 1989) or the 
North Pacific Acoustic Laboratory (NPAL) (Worcester and Spindel, 1997). ATOC seeks 
to measure global wanning related trends based on average large spatial-scale changes in 
the ocean’s temperature distribution. NPAL is a recently proposed semi-permanent 
acoustic laboratory consisting of several vertical and/or horizontal line arrays to be used 
for continuous ocean monitoring and a variety of other experiments. Further introductory 
discussion of issues related to ocean variability and these transmission systems in relation 
to this dissertation are presented in Section B. 

Valid numerical predictions of acoustic travel time variability due to mesoscale 
processes rely critically on realistic environmental inputs. While ideally these inputs 
might come from in-situ measurements, in general the ability to obtain a sufficient number 
of measured synoptic environmental samples to support long range acoustic propagation 
predictions is extremely limited. An alternate approach relies on state-of-the-art ocean 
models to provide the environment input data, and is the approach taken in this disser- 
tation. This is discussed further in Section C. Finally, a summary of the dissertation 
contents is presented in Section D. 
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A. 



HISTORICAL BACKGROUND 



In developing the concept of ocean acoustic tomography within the framework of 
ray theory, Munk and Wunsch (1979) relied on several simplifying assumptions. Two of 
their fundamental assumptions were 

1) travel time changes between widely separated sources and receivers are 
related linearly to the magnitude of the intervening mesoscale 
perturbations, and 

2) ray trajectories will remain stationary upon passage through a mesoscale 
perturbation; only the travel times will vary. 

Mercer and Booker (1983) examined these assumptions and suggested that linearization 
is invalid for dealing with the problem under realistic conditions, particularly when the 
acoustic range exceeds a few hundred kilometers. Their results, based on ray tracing 
through a simulated warm midocean eddy and a cold core Gulf Stream ring (both circular 
with 200 km diameters), suggested that the changes in travel time are nonlinear as a 
function of sound speed. This non-linearity was found to especially affect the late- 
arriving (axial) rays. They further concluded that the assumption of stationary ray paths 
was incorrect, resulting in a significant error compared with the 10-ms precision that was 
expected for tomographic arrays. Overall, their results suggested that the range-dependent 
effects so badly bias the inversions for range-averaged oceans that tomographic results 
might be meaningless. 

Other investigators (Spiesberger and Worcester, 1983) developed comparisons of 
exact and approximate perturbations in travel time and ray geometry due to mesoscale 
perturbations. The three approximations to the exact travel time they considered were 

1) range-dependent perturbation and the unperturbed path, 

2) range-average of the sound-speed perturbation and the unperturbed path, and 

3) range average of the sound-speed perturbation and its associated (perturbed) 

ray path. 

From their study, they concluded that travel-time errors incurred by using either linearized 
or range-averaged travel-time perturbation calculations are small, but not negligible (at 
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most 30% of the exact total travel time). They further concluded that range-dependent 
raytrace calculations would be necessary to achieve accurate inversions. 

By the mid- 1980’s, investigations began to focus on examining the nonlinear bias 
observed in travel time perturbations with a goal of developing systematic methods in 
tomographic techniques to correct for this effect. The nonlinear bias results from the 
differences in ray path geometry (which is a function of the sound speed field) between the 
actual and modeled sound speed fields. Spiesberger (1985a, 1985b) approached the 
problem using numerical prediction within the framework of ray theory, whereas Munk 
and Wunsch (1985, 1987) took an analytical approach, again based on ray theory. General 
results from Spiesberger ’s work suggested that for the non-linear term in a range- 
dependent environment 

1) the non-linear contribution is approximately proportional to the horizontal 
range traveled by the ray, 

2) it is roughly proportional to the square of the magnitude of the eddy 
anomaly, 

3) its value is insensitive to the position of the (isolated) eddies, 
and more generally that 

4) flat, axial rays are less linear than steep rays, 

5) the magnitude and sign of the bias depends on the climatological profile 
used to model the sound speed field, 

6) linear iteration tomographic methods may not be sufficient in all cases to 
correct for the bias, and 

7) the nonlinear bias is small compared with seasonal fluctuations of the travel 
time. 

Munk and Wunsch also found that linear inversions of the range-averaged sound 
speed profile are subject to a significant error arising from a quadratic bias. However, in 
their approach, they found that this quadratic bias is such that the travel time through a 
structured ocean is generally (not always) less than through a smooth ocean with the same 
mean properties. This is referred to as a negative (warm) bias. They also concluded that 
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the bias is generally constant. Thus, what matters for a problem such as detecting global 
warming is not that a warm bias exists, but rather any secular changes in the warm bias. 

Ocean acoustic tomography, when first demonstrated, was performed using 
sources centered at frequencies about 224 Hz and with source-receiver ranges about 300 
km. These conditions were suitable for using a ray theory approach. However, over time, 
there was a growing interest towards refining tomographic techniques for use at longer 
ranges and lower frequencies. Accordingly, to meet the challenges of lower frequency 
and longer ranges, further deep-water tomographic techniques were formed based on a 
combination of ray and mode theory (Munk and Wunsch, 1983), adiabatic mode theory 
(Shang, 1989), a wave-theoretic approach (Athanassoulis and Skarsoulis, 1995), and non- 
adiabatic parabolic equation modeling (Shang and Wang, 1993; Shang, 1997). In general, 
findings with regard to the travel time perturbations are similar to those reported by inves- 
tigators implementing ray theory, with the exception that Shang (1997) reported a predom- 
inant “cold” bias based on full-wave parabolic equation (PE) modeling vice the “warm” 
bias found using analytic ray theory (Munk and Wunsch, 1987). This cold bias is also 
noted in both ray-tracing and PE numerical predictions by Tappert and Tang (unpublished 
1995), in PE numerical predictions by Smith (1992) and Shang and Wang (1996), and PE 
numerical predictions done in support of this dissertation previously presented by Smith 
et. al. (1996). 

Long-range propagation has also been suggested to allow measurements of 
range-averaged sound-speed (and corresponding temperature) profiles, effectively 
filtering out the mesoscale variability in order to observe trends in the mean ocean temper- 
ature (Munk and Forbes, 1989). In contrast to acoustic tomography, where the mesoscale 
variability is the “signal” providing information about the range-dependent sound speed 
mapping, the mesoscale variability is noise in the global measurement problem. If 
constant, the mesoscale bias is not a significant issue as investigators are interested in 
long-term temporal changes in mean ocean temperature. However, if the mesoscale bias 
variability is sufficiently great, then it becomes difficult to separate out the effects of 
global climatic changes from the mesoscale bias. 
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In related research outside the context of tomography, Smith et. al. (1992a, 1992b) 
considered the effects of ocean mesoscale perturbations on long-range acoustic ray 
propagation. It was shown that acoustic rays, particularly those traveling near the axis 
(corresponding to the lowest order modes), exhibited chaotic behavior, i.e., extreme sensi- 
tivity to small-scale perturbations in the environment and/or initial conditions. This again 
suggested that the arrival pattern of the lowest order modes becomes highly complicated 
at long ranges. Thus, chaotic ray motion limits the ability to make deterministic predic- 
tions of long-range arrival time patterns using ray theory. A deep water tomography 
experiment known as SLICE89 (Howe et al, 1991) measured the arrival structure of a 
pulse at a range of 1000 km. The analysis of the resulting data was consistent with the 
aforementioned complexity of the lowest order modes. Further analysis (Colosi et. al., 
1994) strongly indicated that internal waves were responsible for this complexity rather 
than mesoscale fluctuations. However, other preliminary work (Smith, 1992) with Rossby 
waves suggested that such mesoscale features could significantly perturb the near-axial 
arrivals in spite of the adiabatic nature of the mode propagation. 

As this review of travel time perturbation studies suggests, there is general 
agreement that a non-linear bias does indeed exist between travel times observed in a 
structured range-dependent ocean as compared to a smooth range-independent ocean with 
similar mean properties. However, there appears to be divergent conclusions with regard 
to the general sign of the bias. In this dissertation, the emphasis is on not only predicting a 
statistical mean for the bias, but on investigating its variability based on parameters such 
as range, local mode number, and mesoscale perturbation strength. In the deep ocean 
environment, perturbations will be shown to produce a variable overall spreading in the 
arrival structure, primarily producing a delay (cold bias) in the axial arrivals. In the more 
complicated region of the California current, the analysis of the bias remains more 
ambiguous. 
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B. 



ACOUSTIC MONITORING SYSTEMS 



Understanding inherent ocean variability is also important when installing 
"permanent" acoustic monitoring systems. As discussed by Chiu et. al. (1994), there are 
three fundamental acoustic issues which must be carefully addressed when considering 
the locations of permanent sound sources and receivers. These issues include acoustic 
reliability, acoustic stability, and geophysical noise. Acoustic reliability involves 
determining the regions/sites which can be ensonified on a year-round basis based on a 
specified source location. Acoustic stability is concerned with understanding which 
portion of the acoustic arrival pattern is relatively invariant over time; estimation of useful 
travel time series requires stable arrivals. Finally, geophysical noise, or the size of the 
travel time change due to mesoscale features, ocean fronts, or seasonal cycles will directly 
influence the amount of measurement time required before a statistically significant trend 
can be observed. Hence, information regarding the variability of travel time perturbations 
becomes critical to large-scale monitoring projects such as ATOC and NPAL. 

In recent discussions regarding possible site locations for establishing NPAL, the 
SOSUS array site at Pt. Sur, California was suggested as a potential site. Currently, the 
Naval Postgraduate School (NPS) has several sound sources within a 600 km radius of the 
Pt. Sur array site which were specifically deployed to study the California Current system. 
Transmissions from these sources are monitored through a bottom hydrophone at the Pt. 
Sur site. Through the use of coupled ocean-acoustic computer modeling, this dissertation 
specifically investigates issues of acoustic stability and geophysical noise related to the 
acoustic path (474.6 km in length) between one of these sound sources and the Pt. Sur 
receiver. Insights from this variability study might then be potentially useful in array 
configuration design for NPAL, refining tomographic techniques for this region, or for 
assessing travel time variability of ATOC transmissions through this region. 
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c. 



OCEAN-ACOUSTIC MODELING 



In studying the effects of mesoscale variability on underwater sound propagation 
characteristics, it would be desirable to obtain in-situ measurements at sea for the sound 
propagation and oceanographic variables. A complete description of environmental 
conditions for acoustic purposes would require simultaneously sampling oceanographic 
variables of temperature, salinity, and pressure, but this is costly, time-consuming, and 
provides just a snapshot on a single day. Alternatively, coupled ocean- acoustic computer 
modeling provides a cost effective means of studying acoustic propagation through ocean- 
ographic features such as fronts and eddies on a range of spatial and temporal scales. 
Ocean models can also be coupled with atmospheric models to predict future states of the 
ocean in the upper layers. This permits spatial and temporal variability to be examined. 
As with any type of computer model, the model physics must still be validated with 
physical measurements. 

As a result of ever-improving numerical methods, faster computers, and growing 
global data sets, ocean numerical models are becoming increasingly reliable and realistic 
(Semtner, 1995). Thus, state-of-the-art global ocean circulation models can be used to 
investigate geophysical time-scale acoustic variability. A demonstration of this was 
provided by Chiu et. al. (1994) when they presented their results of a ray variability 
simulation of sound transmission from Heard Island to the west coast of the United States. 
A similar ray variability simulation was also provided by Staten et. al. (1996) for an 
acoustic path between Hawaii and Monterey, California. Other coupled ocean-acoustic 
propagation models have been demonstrated by Robinson et. al. (1994) and Newhall et. 
al. (1990), among others. In this dissertation, a global ocean circulation model is coupled 
with a full-wave parabolic equation (PE) acoustic propagation model to study the travel 
time variability in a region approximately 500 km off the California coast The code 
written to specifically interface these specific ocean and acoustic models is a unique 
contribution of this dissertation. 

As discussed by Chiu et. al. (1994), travel time fluctuations yield information 
about the mesoscale processes along the acoustic path. These fluctuations are then useful 
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“signals” for the calibration of eddy-resolving global ocean circulation models. The 
statistics of predicted travel time variability as determined from the coupled ocean- 
acoustic model output and provided in this dissertation may then be checked against 
measured travel time variability for consistency. 

D. DISSERTATION SUMMARY 

The objective of this work is to examine (via computer simulation) the variability 
of modal acoustic travel times due to mesoscale ocean structure and to analyze these 
effects in terms of range, local modes, perturbation strength and travel time variability. 
These results have potential applicability to tomographic and long range acoustic 
thermometry techniques. 

Chapter II, Theoretical Background, reviews the theory relevant to the material 
presented in the subsequent chapters. Overviews are provided of the parabolic equation 
approximation used for acoustic propagation modeling, the PE broadband travel time 
computation method, and the modal travel time computation method Modal Spectrum of 
the PE Field (MOSPEF) (Shang and Wang, 1993). 

Chapter ID, Mode Functions for the Wide Angle Approximation to the 
Parabolic Equation, examines the effect of various parabolic equation approximations on 
the normal mode decomposition of the computed pressure field. Local normal mode 
decompositions were required to analyze modal travel time variability as presented in this 
dissertation. However, in the course of doing this analysis, computations for the modal 
amplitudes of a reference range-independent PE field using standard normal modes 
yielded fluctuating vice constant modal amplitudes for modes higher than about mode 30. 
This led to the development and investigation of approximate “wide angle PE mode 
functions”. These mode functions form a different basis set for modal expansion from that 
obtained using standard normal mode theory. Chapter HI presents numerical implemen- 
tation schemes for approximating these wide angle mode functions, and results comparing 
range- independent wide angle PE field modal decompositions performed using wide angle 
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vs. standard mode functions. While this work was originally pursued as a “side issue” to 
this dissertation, it is perhaps one of the more significant aspects of this research. 

Chapter IV, Deep Ocean Environment, investigates the effects of randomly 
phased Rossby wave baroclinic modes on acoustic travel time and arrival structure in a 
deep ocean environment. Numerical predictions of the acoustic field are obtained from a 
broadband, range-dependent parabolic equation (PE) model using a modeled deep ocean 
environment composed of a superposition of randomly phased Rossby waves to provide 
realistic mesoscale sound speed perturbations. The variability in arrival structure and the 
bias resulting from approximating the range-dependent environment with a range- 
averaged environment is examined for ranges of 500 and 1000 kilometers. 

Chapter V, Coupled Ocean-Acoustic Modeling through a region of the* 
California Current, investigates the effects of a fluctuating frontal region created by the 
California Current on modal acoustic travel time variability. For this phase of research, 
the state-of-the-art global Parallel Ocean Circulation Model (POCM) (Semtner and 
Chervin, 1992), also known as the Semtner and Chervin model, was used to obtain 
realistic input environmental data spanning a one-year time period for a 474 km acoustic 
path. The temporal, spatial and seasonal variability in the individual modal arrivals of the 
first thirty modes is assessed. The mesoscale bias variability is also examined by 
comparing the full-field peak arrival times for the range-averaged environment to that of 
the range-dependent environment. 

Chapter VI, Conclusions, presents a summation of conclusions and identifies 
areas requiring further research. 
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II. THEORETICAL BACKGROUND 



In this chapter, basic background theory relevant to subsequent chapters of the 
dissertation is presented. Specifically, overviews are provided of the parabolic equation 
approximation used for acoustic propagation modeling, the PE broadband travel time 
computation method, and a modal travel time computation method Modal Spectrum of the 
PE Field (MOSPEF) (Shang and Wang, 1993). 

A. PARABOLIC EQUATION APPROXIMATION 

1. Background 

In a 1974 lecture, Tappert (1977, paper based on lecture) introduced the parabolic 
approximation to the wave equation to the underwater acoustics community. This approx- 
imation is valid to much lower frequencies (>2 Hz) than the geometrical acoustic (small 
wavelength) approximation as it retains all the diffraction effects associated with the 
ocean sound channel (Tappert, 1977). In contrast to separation of variables methods, 
which are based on the approximation that the ocean is exactly horizontally stratified, the 
parabolic approximation retains full coupling between local waveguide modes, thereby 
making it valid for more realistic, non-stratified oceans. These are just a few of the 
reasons why use of the parabolic approximation for acoustic propagation predictions has 
become commonplace. 

In Tappert’s original paper, several approximations to the so-called “square root 
operator” of the parabolic approximation were introduced. One of these defines what is 
known as the standard parabolic equation (SPE). Over the last twenty years, there have 
been additional higher order approximations made to the square root operator of the 
parabolic equation. One such approximation leads to the wide-angle parabolic equation 
(WAPE) introduced by Thomson and Chapman (1983). This WAPE approximation is 
based on an operator splitting by Feit and Fleck (1978), and extends the SPE 
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approximation to accommodate wider angles of propagation. It is also generally less 
sensitive to the choice of a reference wavenumber, k Q . Furthermore, error analysis 

indicated that the effects of phase errors are greatly reduced with this operator. 
Subsequently, the WAPE approximation is favored over the SPE for predicting acoustic 
propagation, and (unless otherwise specified) is the approximation used for all acoustic 
propagation computations made in this dissertation. 

2. Theory 



To develop the parabolic equation approximation, first consider a time-harmonic 
acoustic field in a cylindrical coordinate system which can be represented by 

P(r,z,<?,t) = p(r,z, (2.1) 



Substituting this into the wave equation in cylindrical coordinates leads to the Helmholtz 
equation 



LL(» 

rdr\ dr) 



+ 



2 

1 dp 

~~2 ~. ~2 

r d( p 



2 

dp 



dz 



2 +k 0 



2 2 

n (r,z,q>)p 



-4nP 0 $(*-3k t ), 



( 2 . 2 ) 



where k n = — is the reference wavenumber, n(r,z,<0) = — is the acoustic 

u c 0 . c(r,z, (p) 

index of refraction, c Q is the reference sound speed, and c (r,z, (p) is the acoustic sound 

speed. Here, all of the features of the environment are represented within c (r,z, <p). We 
shall neglect density variations in this derivation, but this could be incorporated without 
any loss of generality. For a point source at coordinates (r = 0, z = z s ) , 

= 2Fr 8(l " Z ' )8W ' (2 ' 3) 

The reference pressure level P Q is defined as the pressure amplitude at a reference 
distance of R Q = 1 m. 
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Treating the ocean as a waveguide with a cylindrical coordinate system, acoustic 
energy primarily propagates outward from an acoustic source in the horizontal direction. 
Therefore, the pressure field can be represented by an outgoing Hankel function slowly 
modulated by an envelope function, that is, 

P ( r > z> <P) = V O, z, cp) (*o r )- (2.4) 



In the far-field the Hankel function can be approximated by the first term of its asymptotic 
expansion (Gradshteyn and Ryzhik, 1994) 




(V) “ 



ink Q r 



1/2 ik 0 r 

e 



, k n r » 1 . 



(2.5) 



Using this expression, the relationship between the acoustic pressure p (r, z, cp) and the 
envelope function \| f{r,z, cp) can be expressed as 



P(r,z, cp) 





y(r, z, cp) e 



ik 0 r 



( 2 . 6 ) 



This is the definition of the so-called “PE field function” \)/ , scaled such that at r = R 0 , 
|\jr| = 1 and \p\ = P Q . Substituting this into the Helmholtz equation, assuming 

azimuthal symmetry and dropping the source term, it can be shown that the far-field 
outgoing wave in layered media satisfies (Tappert, 1977; Smith and Tappert, 1993) 

= ( 2 - 7 > 



where the operators 



Q, p = (1 + n + e) 1/2 , 



1 3 ^ 
M- = ~ 2 — » 
k o dz 



e = n -1. 



( 2 . 8 ) 

(2.9) 



( 2 . 10 ) 
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For range-independent environments, it should be noted that Eq. (2.7) is exact and defines 
the solution of the field propagating outward from the source. 

The second order differential equation is now reduced to a first order one, thereby 
allowing solutions to be approximated via a one-way (no backscatter) non-iterative 
marching algorithm 

\\f (r + Ar,z) = exp{ik 0 Ar(Q op - l)}\|/(r,z) . (2.11) 

However, use of the full-band, exact square-root operator is not numerically efficient. 
Subsequently, approximations to Q have been introduced to reduce the computational 

burden. Note also that as approximations to Q are used, the pressure field computed 

from an approximate \|/ substituted into the exact relation Eq. (2.6) also becomes 
approximate. The resulting approximate pressure field is therefore not exactly equivalent 
to the acoustic pressure field satisfying the acoustic wave equation. The significance of 
this will become apparent as the SPE and WAPE mode functions are compared in Chapt. 

m. 

The SPE operator due to Tappert (1977) is obtained using the binomial approxi- 
mation 

Q op = C 1 +V + Z) 1/2 “ 1 + + for |jj.| « 1 , |e| « 1 . (2.12) 

This approximation assumes that variation in the refractive index is small to some degree 
(|e| « 1 ), and is limited to small propagation angles (|fi| « 1 ). The literature suggests 
accurate solutions are limited to a half-beamwidth of 10° - 20° (see e.g., Jensen et. al., 
1994, p.346, or Chin-Bing et. al., 1993, p.62) for the propagation angle. Substituting Eq. 
(2.12) into Eq. (2.7) yields the well-known standard parabolic equation for underwater 
acoustics 




14 



This equation is valid in the far-field, and neglects azimuthal coupling and backscatter. 
Adopting the notation 

Qspe = (^ + 2^ + 2 £ )’ (2-14) 



Eq. (2.13) becomes 

57 = i*o(e sre -l)V. ( 2 . 15 ) 

A higher order approximation to Q introduced by Thomson and Chapman 

(1983) is based on an operator splitting by Feit and Fleck (1978). Commonly referred to 
as the “wide-angle” approximation, it has the form 

Qwape= ( 1 + 1-0 1/2 + [ ( 1 1 e ) 1/2 - 0 • (2-16) 

This operator extended the SPE operator to accommodate wider angles of propagation and 
was also found to be less sensitive to the choice of k Q . The literature suggests accurate 

solutions are limited to a half-beamwidth of -40° (see e.g., Jensen et. al., 1994, p.350) for 
the WAPE propagation angle. However, benchmark testing using the WAPE 
approximation demonstrated the capability to accurately propagate fields for select 
environments with half-beamwidths of greater than 70° (Chin-Bing et. al., 1993). Error 
analysis indicated that for typical deep ocean conditions, the effects of phase errors were 
greatly reduced relative to the SPE operator. For these reasons, use of the WAPE is 
commonplace for predicting sound propagation, and is the approximation used in the 
Navy standard PE model (Holmes and Gainey, 1991). 

These are but two widely-used approximations to the exact square-root operator 
Q defined by Eq. (2.8). Other approximations exist, and use of any of these approxi- 
mations to the exact square-root operator will give rise to some phase error beyond some 
range. In the next section, we shall examine how the wide angle approximation to the 
square-root operator is implemented for numerical predictions in the University of Miami 
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Parabolic Equation Model (UMPE) (Smith and Tappert, 1993), the model used for this 
dissertation work. 

3. Implementation of the PE approximation 

There are three common methods of computing PE solutions: (1) the split-step 
Fourier (PE/SSF) method (Hardin and Tappert, 1973), (2) the implicit finite difference 
(IFD-PE) method (e.g., Lee and Botseas, 1982), and (3) the finite element (FEPE) method 
(e.g., Collins, 1988). UMPE uses the PE/SSF method and the general algorithm behind 
this method is described below. For a more formal development, the reader is referred to 
Smith and Tappert (1993). 

To initiate the computation, a regular spatial grid (Ar, A z) is established onto 
which the environmental parameters (n, p ) are mapped. Next, the marching algorithm, 
Eq. (2.11), requires that the initial conditions for the PE field function \ \f (r Q , z) as well as 

boundary conditions at the sea surface (z=0) and at the maximum computational depth 
(z=z m ax) be specified. The free surface is typically treated as a pressure release boundary, 

thus requiring \(/ (r, z = 0) = 0 . The bottom boundary condition at z - z max must 
approximate the radiation condition (z) —> 0 as z — » . However, since the computa- 

tional depth is finite, the field amplitude is forced to zero in an artificial absorption layer 
which extends well below the bottom depth of interest. Because the PE/SSF algorithm 
numerically solves the one-way wave equation, Eq. (2.7), over all depths, there is no need 
to incorporate a numerical scheme to satisfy the appropriate boundary conditions at the 
water/sediment interface. These boundary conditions are automatically satisfied within 
the PE/SSF algorithm by simply including the changes in the environment at the water/ 
sediment interface within the definition of the acoustic index of refraction, n(r,z). 

Before considering implementation of the marching algorithm further, note that 
the square root operator Q op , Eq. (2.8), is written in terms of two other operators. One 

operator, e , Eq. (2. 10), is simply a multiplication operator in z-space and therefore forms a 
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diagonal matrix on the computational depth grid. The other operator, |i , Eq. (2.9), is a 
function of the second derivative in z-space, and thus forms a tridiagonal banded matrix 
on the computational depth grid. Correspondingly, when implementing the marching 
algortihm, the eigenfunction value in the z-domain at a given depth grid point is coupled 
to values at other depth points. However, since there is a unique eigenfunction for each 
eigenvalue, the corresponding operator in vertical wavenumber space is diagonal. This 
makes it desirable to separate the application of each operator to the domain in which it 
forms a diagonal matrix. The PE/SSF algorithm specifically assumes that these 
component operators, which are in some representation diagonal, may be applied indepen- 
dently. Using the simple relation 

Qop = '-(T op + U op ), (2.17) 

the marching algorithm, Eq. (2.11), may be rewritten as 

V (r + A r, z) = exp {-ik Q ArT op } exp {-ik 0 ArU op } \|/ (r, z) (2.18) 



where T gp and U op come from separating terms derived from |n and £ , respectively, in 

the desired square root operator approximation. This then leads to the following “split- 
step” implementation. 

Once the boundary conditions and initial conditions for (r Q , z) are specified, a 
transformation is made to the fc z -domain followed by a multiplication of the & z -space 



operator e lk ° Ar ^ op ^ . The result is then transformed back to the z-domain and is 

- ik Q &rU (z, r) 

followed by a multiplication of the z-space operator e . The final result is the 

field function at r + Ar. Therefore, the approximate solution for T is marched out in 
range according to the relation 



'E (z, /" + Ar) 



-ik 0 ArU op (z, r) 
e 



xIDFT{e 



-ik 0 &rt op (fe z ) 



x [DFTi'V (z, r))] } 



(2.19) 



where for the wide angle PE approximation 
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op ^ z } 







(k 1 


2“ 


= 1 - 


1 - 










K k 0y 





1/2 



and 

U op^ r ^ = -[«(^ r ) - !] • 

UMPE utilizes a FFT algorithm to implement the DFT in Eq. (2. 19). 



B. PE BROADBAND TRAVEL TIME COMPUTATION 



( 2 . 20 ) 



( 2 . 21 ) 



The calculation of the time domain arrival structure is straightforward using the PE 
code. For a broadband acoustic field, the time-harmonic acoustic field initially assumed in 
Eq. (2.1) represents a single (CW) component of the general field. Separate PE calcula- 
tions are then made to compute acoustic fields for each desired frequency over the 
bandwidth of interest. Using discrete time, / = nAt, and discrete frequencies, / = lAf, a 

total of N frequencies in the bandwidth of interest where A/A/ = ^ , assuming azimuthal 

symmetry and adopting the notation 

P(r,z,nAt ) = p l (r,z)e , (2.22) 



the representation of the field in the time domain can then be represented as 



N- 1 



P (r,z, nAt) = IDFT[X(lAf) Pl (r,z)] = ^ X (lAf) P[ (r, z) e 






(2.23) 



1 = 0 



where X (/A/) represents the complex source spectrum. For all work presented here, a 
Hanning window is used to represent the source amplitudes providing a simple 
representation of a time-domain pulse. 

In UMPE, arrival time structures are computed at user-designated fixed ranges, i.e. 
at r = R , using the value of the frequency-dependent range-scaled PE field function 
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Consequently, the time domain is heterodyned around the value t Q = R/c Q . Accord- 
ingly, arrival times are given in terms of “reduced time” or t - t Q . To reduce the computa- 
tional burden, the UMPE model also heterodynes the signal by shifting the center 
frequency f c to d.c. Adopting the notation/, = /A/, f c = l c Af, f = l u Af and 

t 0 = n Q At , where f l ,f c , and f u represent the lower, center and upper frequency for the 
bandwidth of interest, respectively, the inverse Fourier transform actually computed by 



yielding P' R ( z , (n - n Q ) At), the set of complex pressure values in depth-(reduced) time 
space for the designated range. 

C. MODAL ARRIVAL TIME COMPUTATION 

To compute individual modal arrival time structures, the PE solution must first be 
decomposed into normal modes. Normal mode decompositions based on the technique of 
separation of variables assume that the ocean is nearly exactly stratified horizontally. In 
general, range and depth dependent spatially continuous modes exist which are unique 
harmonics of the ocean waveguide. However, for cases where the ocean acoustic 
waveguide properties vary with range, these modes are not simply functions of depth 
where separation of variables can be used. Subsequently, normal mode theory has been 
generalized for that of an “almost stratified” medium. In ocean acoustics, this theory was 
first set forth by Pierce (1965) and Milder (1969). It has become common practice for 
the modes of an exactly stratified ocean to be used as an orthonormal basis set for any 
general ocean environment. In ocean acoustics, this orthonormal basis set is referred to as 



UMPE is 




(2.24) 
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the set of “local normal modes”, and computed as a function of range using the local 
environmental properties. The result of this approach when applied to the wave equation 
for a range-dependent environment is a set of coupled differential equations which govern 
the range variations of the modal amplitudes. 

Consequently, it is seen that in a range-dependent environment, the local normal 
modes do not propagate independently of each other but rather interact continuously with 
neighboring modes. If, however, the range-dependent changes in sound-speed profile are 
small and gradual, an adiabatic approximation can be assumed. That is, coupling between 
modes may be neglected and each normal mode is assumed to propagate independently 
while it progressively adapts to changes in the environment. Under adiabatic conditions, 
the largest proportion of wave energy will remain in the same-numbered mode as it propa- 
gates through the changing environment. However, it has been shown by Desaubies et. al. 
(1986) that the degree to which the adiabatic approximation can be used is very sensitive 
to the shape of the sound-speed profile, frequency of propagation, and mode number. 
Their work points out that when coupling cannot be neglected, there is no simple 
expression for the modal travel times. 

In an attempt to address this difficulty, Shang and Wang (1993) proposed a method 
of computing modal travel time for range-dependent, strong mode coupling environments 
which is known as the modal spectrum of the PE field (MOSPEF). It is this method 
which is used to compute the modal travel times presented in Chapt. V of this dissertation. 
In MOSPEF, coupled modal amplitude equations are avoided by first computing the total 
acoustic pressure field solution using the PE method, such as that described in previous 
sections. Next, using local normal mode functions 4 / m (z;r, co) computed from a normal 

mode code such as KRAKEN (Porter, 1991), range-dependent modal amplitudes and 
phase are then calculated. 

Using the PE solution for the two-dimensional range-dependent acoustic field, this 
field is expressed as a summation of local normal modes 

oo 

Pp E (r,z;C 0 ) = 2 ( r ;®) ( z;r > ®) ( 2 - 25 ) 

m = 1 
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where 4 / m (z) is the local mode solution of the range-dependent Sturm-Liouville 
eigenvalue equation 

« + W" fe - <] (*) = 0 (2.26) 

dz 



evaluated at the local range r and subject to appropriate boundary conditions. In this 
equation, represents the separation “constant” (which locally assumes a range- 
independent environment) and 'V m (z) denotes the specific mode function associated with 
this separation constant. 

Taking advantage of the orthogonality of the mode functions leads to 

( r ; c °) = jp PE (r,z-,0i) y i , * m (z;r,(0)dz, (2.27) 



or, by substituting Eq. (2.6) to put Eq. (2.27) in terms of the PE field function 

3> m (r;co) = \P 0 ^^ PE {r,z\(Si)^* m {z\r,^) d^e k ° r (2.28) 

- 1 / 2 . . . «V - 1 / 2 ,. . . , 

= r A m (r;co) e = r |A m (r;co)|e 

The complex modal amplitude is given by 

A m (r;(a) = jP 0/ JF 0 \\r pE (r, z‘, (£>)¥* m (z;r,(£>) dz (2.29) 



and the modal phase is given by 

<P m O;a>) = k 0 r + Arg[A m (r;a)] . 



(2.30) 



Shang and Wang refer to the quantity | A m (r;co) | as the modal repopulation character. 

The modal arrival times can be computed similarly to the broadband computation. 
Using discrete frequencies with the source spectrum denoted by X (/A/), the acoustic 
pressure field is found from the PE pressure field solution by substituting Eqs. (2.25) and 
(2.28) into (2.23), thus yielding 
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Hr, z, nAt) = i X Z x V*f> r ~ 1/2 \A m (r;lAf) \V m (z;r, /A/) e 



Af-l/ 




y 



or in terms of the modal pulse 




1 = 0 



If the mode functions are approximately invariant over the frequency bandwidth of 
consideration, such that 



where (0 [ , co c , and (£> u are the lower, center, and upper frequencies of the bandwidth under 
consideration, then Eq. (2.32) can be simplified to 



This simplification reduces the computational burden of computing an FFT for each depth 
grid point. 

In the dissertation work presented here, modal arrival time structures are computed 
only at designated fixed ranges. The complex modal amplitude A m (R; go) was computed 
from Eq. (2.29) using the individual frequency-dependent mode functions (rather than the 
approximate mode function 'i' m ( z,R , (0 C ) when Eq. (2.33) is valid) as computed by the 

KRAKEN normal mode code. When computing the modal pulse, Eq. (2.32) or Eq. (2.34), 
the approximate center frequency mode function was used only if Eq. (2.33) was valid. 
The signal is heterodyned about the center frequency and the overall phase factor is 
neglected, giving results in terms of “reduced time”. The source spectrum used by UMPE 
is a Hanning amplitude window. Thus, the Fourier transform actually computed was 




(2.33) 




1 = 0 



p 'R,m(z,(.n-n Q )At) = 




(2.35) 
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when Eq. (2.33) was valid or 



P'R, m (z,(n-n 0 )At)=jj X C z;R , /A/) )A m (i?;/A/) « 






/=/, 






(2.36) 



otherwise. The final result yields P'/{ ( w (z, (n - n Q ) At), the set of modal complex 
pressure values in depth-(reduced) time space for the designated range. 
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III. MODE FUNCTIONS FOR THE WIDE ANGLE APPROXIMATION TO THE 

PARABOLIC EQUATION 



With a goal of further developing modal ocean acoustic tomography, Shang and 
Wang (1993) proposed using the modal spectrum of the PE field for computing modal 
travel times under strong mode coupling environments. Instead of dealing with the 
coupled equations governing the complex modal amplitudes, they obtained the coupled 
modal amplitudes and phases from the spectral expansion of the full wave field given by 
the PE solution with local normal modes as the basis set. Their results, obtained using a 
2-D wide angle implicit finite-difference (IFD) PE code (Lee and Botseas, 1982) and 
normal mode code KRAKEN (Porter, 1991), presented modal travel time perturbations for 
the first ten modes due to a single mesoscale eddy. 

In a similar type approach, modal decompositions of wide angle PE fields were 
computed in support of the investigations presented in later chapters. However, computa- 
tions for the modal amplitudes of a reference range-independent wide angle PE field using 
standard normal modes erroneously yielded fluctuating rather than constant modal ampli- 
tudes for modes higher than about mode 30. This led to the new development and investi- 
gation of “wide angle PE mode functions” presented here. The results of Shang and Wang 
do not contradict these results. Rather, their analysis was limited to only a small number 
of low-order modes which do not suffer significantly from this mismatch. It should be 
noted, however, that their analysis did not employ the correct modal basis set, as will be 
shown in this chapter. 

In that PE-based methods are full wave methods providing the total acoustic field 
on a computational grid, they do not readily supply information regarding the propagation 
of individual modes. This modal information might assist in interpreting acoustic propa- 
gation behavior in a complex environment, or be useful for determining the relative 
adiabatic nature of a particular environment. Thomson (1993) addressed this issue for the 
standard parabolic approximation. In his work, he described a PE-based spectral method 
suitable for modal analysis. Specifically, he demonstrated that the field p in a waveguide 
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satisfying the acoustic wave equation is exactly related to the field V| / satisfying the SPE 
approximation. For this case, the modal amplitudes and wavenumbers of \jr can be deter- 
mined exactly, and the corresponding amplitudes and wavenumbers of the normal modes 
of the acoustic field p are then obtained by simple mapping rules. 

For other higher order parabolic approximations, such as the WAPE and Pade' 
expansions, the modal eigenfunctions and associated eigenvalues of \jr are distinct from 
those of the acoustic wave equation. In that the WAPE is a common implementation for 
PE models using the split-step Fourier technique and is expectecf to yield more accurate 
results than the SPE, a better understanding of its “normal modes” is desired. 

In this chapter, the parabolic approximation to the wave equation is examined 
within the context of normal mode theory. Building on the PE theory presented in Chapt. 
II, approximate mode functions for the split-step Fourier (SSF) WAPE algorithm are 
developed. These mode functions are then used to decompose range-independent sound- 
pressure fields computed using the WAPE approximation. The resulting modal coeffi- 
cients and eigenfunctions obtained using the WAPE mode functions are compared with 
those obtained using standard normal mode theory. It is important to note in what follows 
that use of the SSF algorithm does not introduce any additional approximations in terms of 
the corresponding modal description. 

A. THEORETICAL DEVELOPMENT 

Normal mode decomposition theory uses the technique of separation of variables 
and is based on the approximation that the ocean is exactly horizontally stratified. 
Considering only modes within the water column, and assuming constant water density, 
azimuthal symmetry, and range independence (sound speed varies only with depth), the 

time-harmonic acoustic field p (r,z) e l(0t at r > 0 due to a point source at z = z and 
r = 0 satisfies the two-dimensional Helmholtz equation 
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(3.1) 



: » V<z)P = °' 



Seeking solutions in the form p (r, z) = O (r) T* (z) yields the well-known 
depth-dependent modal equation 

-^ m <z) + [*o 2 "^)-CK« =o 

dz 



(3.2) 



which is subject to appropriate boundary conditions. In this equation, K represents the 

separation constant and 'F (z) denotes the specific mode function associated with this 

separation constant. The complete sum of normalized modes is necessary to represent an 
arbitrary pressure field. Accordingly, the pressure is expressed as 

co 

P(r,z) = X 4>„(r)'F m (z). (3.3) 

m = 1 

Returning to the relationship between the acoustic field and the envelope function, 
Eq. (2.6), an alternate approach employing the separation of variables technique to solve 
the standard parabolic equation, Eq. (2.15), yields an expansion into PE modes. Hereafter, 
these particular modes will be referred to as SPE modes. In this case, the range- 
independent expansion is shown by Desanto (1977) to be 

oo 

Y( r >z) = X a j u j( z ) exp 0‘y) (3.4) 

7=1 

where a. are the SPE mode amplitudes, u. (z) are the SPE mode functions, and Sj are the 
SPE modal wavenumbers. 

Thomson (1993) has shown that the depth-dependent SPE mode functions Uj (z) 

satisfy an eigenvalue equation comparable to that for the normal mode functions of the 
acoustic wave equation, Eq. (3.2) 
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(3.5) 



-^2 u j ( z ) + [V" 2 ^ ~ rf] U j ^ = °> 

dz 



where the SPE modal wavenumbers Sj are given by 



s i 



rf -* 0 

2k n 



(3.6) 



and the horizontal wavenumbers K- are the corresponding normal mode eigenvalues of 

Eq. (3.2). The important issue here is that the mode functions for the SPE approximation 
exactly match those which solve the acoustic wave equation. Accordingly, these mode 
functions may be found using a normal mode code such as KRAKEN (Porter, 1991). As 
we shall see shortly, this is not the case when using the WAPE approximation. 

We now seek approximate solutions to Eq. (2.7) replacing Q by the approx- 
imate propagator function Q WAPE and assuming 

oo 

V(r,z) = £ r „(r) v „(z). (3.7) 

n = 1 



Here, T n ( r ) represents the WAPE range dependent modal amplitude and v n (z) is the 
WAPE mode function. This expansion for \\i(r,z) leads to the modal equation 



v„U)-F-rr4-r„(r) = 0. (3.8) 



v„(z) 



ITi i 


1/2 -] 


+ 

1 

1 


l *o dz) 





r Ar)dr " 



Setting the separation constant equal to z‘P n , it can be shown that T n (r) has the simple 
range dependence 

r n( r ) = & „ ex P( f P n r ) ( 3 - 9 ) 



where b n is the modal amplitude and (3 n is the associated WAPE modal wavenumber. 
Using this separation constant, the WAPE depth-dependent modal equation becomes 
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1 + 



1 d 

,2 , 2 

k Q dz 



2 \ 



1/2 



+ (n(z) - 1) - 






1+ r 

K 0j 






V„U) =0. 



(3.10) 



To develop approximations for the WAPE modes, the square root operator is 
approximated using first order and second order Taylor series expansions. This leads to 
the first order approximate depth-dependent WAPE modal equation 

,2 -r (3 '• 



u> + 

dz 



n(z) ~T 

K 0 J 



v„(z) =0, 



(3.11) 



and to second order 



J__d_ 

l 8 k^dz 



Id 2 . . P/ 

o — 7 - n (z) ~T~ 

2k 0 dz k o 



v„(z)=0. 



(3.12) 



where P w ' = + k Q . Note that the terms “first order” and “second order” refer to the 

relative order of the Taylor series employed. It does not relate to the order of accuracy 
used during numerical implementation. The first order approximation is essentially a 
return to the small angle approximation, and is not expected to provide highly satisfactory 
results in numerical implementation. However, it is included in the analysis to gain 
insight into the nature of the WAPE mode function and its relationship to the standard 
normal mode function. 

The eigenfunctions satisfying Eqs. (3.10) - (3.12) no longer match those normal 
mode eigenfunctions associated with the acoustic wave equation. The orthogonality (i.e. 
uniqueness) of the WAPE modes can be examined as follows. Here we shall only 
consider the first order approximation, although the orthogonality of the WAPE modes can 
be shown to be more general. Multiply the first order approximate depth-dependent 
WAPE modal equation, Eq. (3.11), for v n (z) by v m (z) and for v m (z) by v n (z) . 

Subtracting these two equations yields 
2 2 

= 2k 0 ( P* “ Pm) V « (*) V m U) ( 3 ‘ 13 ) 

dz dz 
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or 



lL V m {z) Tz Vn{z) -V " (z) ^ V - (z) ) = 2fc 0 (P n -PJ v „(2) v m U) • (3.14) 



Integrating over z and applying the boundary conditions for an ideal waveguide (pressure 
release surface and pressure release bottom), Eq. (3.14) reduces to 



The assumption of a pressure release bottom boundary condition in this case is considered 
appropriate as a penetrable seafloor at low grazing angles reflects similarly to a free 
surface (Jensen et. al., 1994, p. 108). Note also that the assumption of a pressure release 
bottom does not greatly alter the resulting modeshapes when considering only deep, 
trapped modes corresponding to refracted-refracted (RR) paths with no bottom 
interactions, as will be the case for the remainder of this development and subsequent 
numerical implementation. Since each modal wavenumber is unique, i.e. 
if m*n. 



Therefore, the first order approximate WAPE modes form an orthogonal basis set. 

To determine the approximate WAPE modal amplitudes, Eqs. (2.4), (2.5), (3.7) 
and (3.9) are used to establish the far-field WAPE modal solution for the computed WAPE 
field p in a layered waveguide 




(3.15) 



,D 

J 0 v n ( z ) v m ( z ) dz. = 0 for m * n. 



(3.16) 




(3.17) 



where (3/ = P„ + V 
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B. 



NUMERICAL IMPLEMENTATION 



Finite-difference techniques are used to numerically approximate the wide-angle 
PE mode functions developed in the previous section. Numerical schemes are developed 
for the first and second order Taylor series approximations given by Eqs. (3.11) and (3. 12). 

For the first order approximation Eq. (3. 1 1), a numerical scheme closely following 
that developed by Jensen et al. (1994) for normal modes to the acoustic wave equation is 
developed. Defining a grid of N equally-spaced points over the depth interval 0 < z< D , 
gives the depth points Zj = jh ,j = 0,1,... jV, where h is the mesh width given by D/N. This 

then allows us to define the following 0^/z 2 J difference approximations for the first and 
second derivatives based on forward, centered, and backward difference approximations, 
respectively. 



•/•V-Wsi-' •■)•-]!' 

’i^ + o(h 2 ). 



V. , - 2v, + V; 

v." = Jzl L 

j , 2 



(3.18) 



(3.19) 






(3.20) 



In the water column, the approximate mode functions must satisfy Eq. (3.11), 
which in terms of the difference approximations becomes 



v . _ j + [ — 2 + h 2&q 



[ 4 > 



-B ' 

r n 



v. + v. + 1 = 0, j = 1, 2, ...N- 1 . (3.21) 



As discussed previously, the top and bottom boundary conditions are assumed to be 
pressure release, giving v Q = 0 and v N = 0 . Collecting these equations together, we 

obtain the eigenvalue problem, 

A((V)v = 0 (3.22) 
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where A is a symmetric tridiagonal matrix of order N - 1 with the diagonal elements defined 
by 




- 2 + h 2k, 



o 



~ (0 
c ( Zj ) 




,j = 1, 2, ...N- 1 , 



and the off-diagonal elements defined by 

e. = 1 ,j = 1,2, ...N-l. 



(3.23) 



(3.24) 



To develop a numerical scheme to solve the second order approximation Eq. 

(3.12), the following o(^/i 4 J centered difference approximations (Godunov and 
Ryabenkii, 1987) are employed 



IV 

v j 



1 r 1 

*, 4 L 6V 



0 13 28 13 n 

3 + 2v /-2 - T V i + —Vj- -v j + l + 2v, 



i * 2 6 V -/ + 3 ]’ 



(3.25) 



„ _ 1 r 1 4 54 i i 

j A 12 v J-2 + 3 V j-i 2 V j 3 V j +l 12 V -' + 2 J' 



(3.26) 



Substituting into the governing equation 

v' V -«o''"-8^(^ r |3„')v = 0 (3.27) 

and using a mesh similar to that defined for Eqs. (3.18)-(3.24) gives the resulting equation 
which must be satisfied in the water column. 

Vs + (- 12 - 2hIk l )v j + ( 39 + 32 h 2 k l) v j-< 

+ (- 56 - 60 h\l + 48/ i 4 fc^( ^ + ( 39 + 32* 2 *j)v. , , 

+ [-l2-2h 2 k 2 o y j+2 + v j+3 = o, j = 3, 4, ...N-3. (3.28) 

This leads to a seven diagonal banded matrix. Again, pressure release boundaries are . 
assumed, giving v Q = 0 and = 0 . For j = 1, 2, N — 2, N - 1 , the pressure release 
assumption leads to the following equations, 
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( 3 . 29 ) 



(- « - 58 h 2 k] + 48 hX{ ?§) - P„')) v, + ( 38 + 32A 2 tJ Jv 2 
+ (- l2-2h 2 k 2 ^jv 1 + v 4 = 0, j = 1, 

( 38 + 32 hhl)v t + (- 56 - 60 hh\ + - p„'))v 2 

+ (39 + 32A 2 ^)v 3 + (- 12-2k 2 kl)v 4 +v 5 = 0, j = 2, (3.30) 

5 + ( ' - 12 ■ - 2 * 2 ^)v ».4 + (39 + 32 A 2 tj)v„. 3 
+ f- 56 - 60 h 2 kl + 48*X( -|L - P„'))v w _ 2 
+ ( 38 + 32A 2 *o)v w _ , = 0, j = N- 2, (3.3 1) 

v„_ 4 + (- 12 - 2A 2 *o) v at- 3 + ( : 38 + 32 h 2 k\)v N _ 2 (3.32) 

+ (-44 - 58f ! 2 ^ + 48/, 4 ^(^-P;))v Ar _, = 0, j = N- 1 

In the following section, results from these numerical approximations will be shown. 

The eigenvectors of the matrices defined by Eqs. (3.23)-(3.24) and (3.28)-(3.32) 
represent the first order and second order Taylor series approximations, respectively, to the 
WAPE mode functions v n (z) in Eq. (3.16). The matrix eigenvalues represent the 

associated wide-angle PE modal wavenumbers p n . 

To determine the modal amplitudes, it is useful to write Eq. (3.17) as 

S r P (r ’ z) “ ( l£k) 5X v n W ex P ( J 'Pn V ) = 2>n (r) V„ (Z) (3.33) 

where A n is the desired modal amplitude. This expression is general and does not depend 
on which approximation to the WAPE is used. Note that we have accounted for 

cylindrical spreading losses by scaling p(r,z ) by Jr. Then assuming normalized modes 
and invoking orthogonality, it can be shown that at a given range the modal amplitude is 
given by 
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(3.34) 




In the following section, modal amplitudes for a range-independent, open ocean 
environment are computed. 

C. NUMERICAL RESULTS 

The numerical schemes developed for the first and second order approximations to 
the WAPE mode function were tested assuming a range-independent, deep ocean 
environment characterized by the Munk canonical sound speed profile (Munk, 1974) 



where 8=0.0057, the scaled depth variable r| =2(z-Z ax j s )/B, and the reference sound speed 
c Q was chosen as 1490 m/s. Axis depth, z^, was assumed to be -1.0 km. The bottom 
depth was assumed to be -4.5 km. 

The acoustic pressure field was computed out to a range of 100 km using the 

University of Miami PE (UMPE) model (Smith and Tappert, 1993) employing the wide 
angle approximation. Input run parameters are summarized by Table 1. As was discussed 
in Chapt. II.A.3, the computational depth includes the water column of interest (bottom 
depth listed) and a “finite-depth fluid bottom halfspace”. This halfspace extends to the 
computational depth and is terminated with an artificial absorption layer to meet the 
radiation condition. The resulting pressure field was then decomposed into standard 
normal modes using KRAKEN (Porter, 1991), and first and second order approximation 
WA modes by finding the eigenvectors of the matrices defined by Eqs. (3.23)-(3.24) and 
(3.28)-(3.32) respectively. 

These mode functions were then used to compute the respective modal amplitude 
coefficients at 1 km intervals out to 100 km using Eq. (3.34). For a range-independent 
environment, these coefficients should remain constant with range. Plots of selected 
modal coefficient amplitudes for the three different mode functions are provided by Figs. 




(3.35) 
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Parameter 


Value 


Reference sound speed, c re j 


1490 m/s 


Frequency, f 


100 Hz 


Source Level 


183 dB re 1 ^Pa 


Source Spectrum 


Hanning amplitude window 


Source Depth 


1000 m 


Bottom Depth 


4500 m 


Computational Depth 


8000 m 


Transform size 


8192 


Depth grid step 


0.9765 m 


Range grid step 


10.0 m 



Table 1 . Acoustic propagation model (UMPE) input run parameters used to compute 
WAPE solution for the modeled range-independent deep ocean environment. 

1-3. Modes were selected so as to minimize overlapping line segments on the plots. Only 
non-bottom interacting modes were considered, as the numerical scheme for the second 
order WAPE mode functions currently only treats the bottom as a pressure release surface. 

For the lower modes depicted by Fig. 1, all three methods provide satisfactory 
demonstration of range-independence. However, the first order WAPE mode function 
approximation starts to show some small fluctuations for mode numbers greater than 15. 
In Fig. 2, the first order approximation clearly breaks down as the mode number increases. 

Additionally, the standard normal mode decomposition of the WAPE computed 
field begins to show noticeable fluctuations for modes above mode 35. Only the second 
order approximation demonstrates satisfactory range-independence. As the mode number 
is increased further, the fluctuations in both the SPE and first order approximations 
continue to increase. This is illustrated in Fig. 3. Again, the second order approximation 
continues to demonstrate the range-independent nature of the input sound-speed 
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Figure 1. Modal coefficient amplitudes as a function of range for modes 1,8, 15, 20 and 
25 of a range-independent Munk canonical sound speed profile. Mode shapes were 
determined using standard normal mode theory and first and second order approximations 
to the wide angle mode function. 
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Amplitude Amplitude Amplitude 



Standard normal modes 





Range (km) 

WA modes 2nd order approximation 




Range (km) 

Mode 30 b Mode 35 x Mode 40 + Mode 45 * Mode 50 o 

Figure 2. Same as Fig. 1 except for modes 30, 35, 40, 45, and 50. 
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Amplitude Amplitude Amplitude 



Standard normal modes 






Range (km) 

Mode 55 b Mode 60 x Mode 70 + Mode 72 * Mode 74 o 

Figure 3. Same as Fig. 1 except for modes 55, 60, 70, 72, and 74. 
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environment. This pattern was observed up through mode 100, the highest mode 
computed in this work. 

The differences are relatively small between the mode shapes computed using 
standard normal mode theory and the second order approximation to the WAPE mode 
function derived in this chapter. Figure 4 illustrates the difference between these two 
mode functions for mode 70. For this particular sound speed profile, the mode shapes 
derived using the two methods tended to converge at low mode numbers, and slowly 
diverge as the mode number increased. That such small differences in mode shapes would 
generate significant errors in the mode amplitude decomposition is both an interesting and 
important finding of this disseration. 

Finally, in that a portion of the work in this dissertation uses mode functions to 
compute modal travel times, a comparison of modal travel times computed using standard 
normal mode functions and the second order approximate WAPE mode functions is 
depicted by Fig. 5. For the results presented in Chapt. V, standard normal mode functions 
were used to compute modal travel times. This choice was made for two primary reasons. 
The first reason is that we were interested in statistical trends in the data rather than 
absolute times. The relatively small difference in travel time results between the two 
methods as represented by Fig. 5 was deemed acceptable for these trends. The second 
reason is that the numerical implementation for second order approximate wide angle 
mode functions developed in this chapter does not yet allow transmission across an 
interface or the incorporation of a density contrast, such as was required for the water 
column / bottom layer interface in Chapt. V. 

D. SUMMARY 

This chapter has examined how the wide angle approximation to the PE square 
root operator alters the associated mode functions in a normal mode decomposition. 
Specifically, first order and second order Taylor series approximations to the WAPE mode 
functions were developed and compared to SPE mode functions. While the field \|/ 
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Figure 4. Comparison of mode 70 as computed using standard normal mode theory and 
the second order approximation to the wide angle mode function. The left hand side 
depicts the modeshape over the entire water column. The right hand side provides an 
enlargement of the 3200-3400 m depth to more clearly depict the difference in mode 
shape. 
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Figure 5. Comparison of arrival times computed using standard normal modes 
(KRAKEN) and second order approximate WAPE modes. Data illustrates the arrival 
times computed for data set 36 (95Jul07) as described in Chapt. V. The x,o plots the 
centroid of the arriving acoustic energy, and the errorbar plots the first standard deviation 
of energy distribution about the centroid. The mean difference in arrival times over the 30 
modes was 2 (±15.9 ) msec., with the positive mean indicating that the WAPE mode 
functions tended to predict a slightly earlier arrival (for this environment) relative to the 
decomposition estimates using standard normal modes. 
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satisfying the SPE approximation is exactly related to the field p which solves the acoustic 
wave equation, this is not the case for the WAPE approximation. 

Numerical schemes were developed to compute the approximate WAPE mode 
functions and wavenumbers. The computed mode functions were then used to decompose 
a typical deep ocean range-independent sound-pressure field obtained using the WAPE 
approximation. The resulting modal amplitudes obtained using the WAPE mode 
functions were compared with those obtained using a standard normal mode 
decomposition of the WAPE field. The second order approximation to the WAPE mode 
function showed noticeably less fluctuation in the computed modal amplitudes for the 
modeled range-independent environment in comparison to the standard mode shapes 
above mode 35. The first order approximation to the WAPE mode function, although 
simpler to implement numerically, did not provide an adequate representation of the 
range-independent nature of the modeled environment. 

The purpose of this chapter was to highlight the need for a different set of mode 
functions in normal mode expansions for parabolic approximations to the acoustic wave 
equation other than the SPE approximation. While the second order Taylor series approx- 
imation to the WAPE mode function derived in this chapter showed promising results for 
the simplified deep ocean environment, much remains to be done to develop proper 
numerical treatment of interfaces, bottom boundaries, and density gradients. For the 
modal arrival analysis accomplished in support of Chapt. V, the KRAKEN normal mode 
theory was used to compute modal arrival times. 
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IV. TRAVEL TIME BIAS DUE TO OPEN OCEAN ROSSBY WAVES 



In this chapter, the effects of randomly phased baroclinic Rossby waves on 
acoustic travel time and arrival patterns are investigated. Rossby waves are characterized 
by a small sea-surface signature (10 cm or less), slow propagation speeds (of order 

10 cm s' 1 or less), and long wavelengths (hundreds to thousands of km). It has been 
speculated that these waves are responsible for the westward intensification of circulation 
gyres (Chelton and Schlax, 1996). It will be shown that the perturbations resulting from 
this random mesoscale structure produces an overall spreading in the arrival structure, 
primarily producing a delay (cold bias) in the axial arrivals. 

While other investigators using various methods of analysis have reported a non- 
linear mesoscale bias which particularly affects the late arrivals (e.g., Mercer and Booker, 
1983, Smith, Brown and Tappert, 1992a, 1992b, Athanassoulis and Skarsoulis, 1995), an 
extensive literature search found no full-wave field studies which specifically address the 
variability in the observed bias for a field composed of randomly phased Rossby waves. 

A further understanding of this bias effect and its variability could provide improved 
predictions of acoustic pulse propagation and provide insights into methods which correct 
for this bias in ocean monitoring efforts such as ATOC. 

A. THE MODELED OCEAN ENVIRONMENT 

Earlier investigations (Smith, Brown, and Tappert, 1992b) with Rossby waves 
suggested that such mesoscale features could significantly perturb the near-axial propa- 
gation. As a follow-on to this previous work in studying the variability in the observed 
perturbations, a similar modeled ocean environment was used in this analysis. Major 
features of the model along with analytic expressions for both the background and pertur- 
bation sound-speed fields are provided below. Expressions for the sound-speed fields 
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assume the Brunt-Vaiasala (buoyancy) frequency has the exponential form 






(4.1) 



where B = 1 km, N 0 = 3 cycles/h, and the depth z increases upward. The ocean surface 
and bottom are assumed to be flat, perfectly reflecting interfaces at depths z=0 and -4.5 
km, respectively. The range-independent background sound-speed field consistent with 
this exponential form of the buoyancy frequency can be defined by the Munk canonical 
profile (Munk, 1974) 



where £=0.0057, the scaled depth variable T| =2(z-z ax i s )/B, and the reference sound speed 
c 0 is chosen as 1490 m/s. Axis depth, z^, is assumed to be -1.0 km. 

To simulate realistic mesoscale ocean features, the modeled environment contains 
a superposition of several modes of randomly phased baroclinic Rossby waves. Chiu and 
Desaubies (1987) have shown that this superposition is appropriate for modeling an open 
ocean environment away from strong boundary currents. The wave structure produces 
vertical water displacements of the water column, C,(x,y,z,t). The coordinate system is 
aligned so that x, y, and z increase to the east, north, and upwards respectively. The sound- 
speed perturbation 5 c(z,x) due to this mesoscale structure can be computed from a 
knowledge of N(z) and the vertical water displacement C,(x,y,z,t) (Munk and Zachariasen, 
1976). Smith et. al. (1992b) have further shown that this sound-speed perturbation can be 
computed from 




(4.2) 



Sc (z, x) 




(4.3) 



with 
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z/ B 

where £ . (z) = C-Be' . In Eq. (4.3), f Q is the Coriolis parameter, which we have 

assumed to be 1 cycle/day corresponding to a latitude of 30° , fi =24.5, g=9.8 m/s 2 is the 
gravitational acceleration, and V 0 is the perturbation strength parameter in m/s. Each 
mode is assigned a weight Aj, wavenumber kj, and random [uniform on (0,2 ti)] phase a ; . 
In Eq. (4.4), the mode parameters Cj are selected to satisfy the bottom boundary condition 

Y o[ C j Be ~ H/B )~ [ F 0 ( CjB) /J Q ( CjB) } J Q [ CjBe~ H/B ) = 0. (4.5) 

In the work reported here, the first four baroclinic modes are included in the modal 
expansion. General experience has shown the lowest few modes contain most of the 
energy (Munk, Worcester, and Wunsch, 1995). Accordingly, mode parameter values were 
assigned so that the normalized modal energies dropped exponentially with increasing 
mode number. Assigned mode parameter values are Xj=2n/kj=400, 350, 300, 250 km and 
A ; = 1.0000, 0.5157, 0.3276, 0.2236 for y=l,2,3,4 respectively. These parameters then 

provide normalized modal energies which scale as ( Aj/kj ) 2 of 1 .00: 0.20: 0.06: 0.02. 

The work of Smith, Brown and Tappert (1992b) reported that the threshold pertur- 
bation strength for the onset of easily detectable chaotic ray behavior appeared to lie 
between V o =0. 1 25 m/s and V^=0.25 m/s, perturbation strengths typical of the mid- to low- 
latitude regions of the worlds ocean. To explore this same perturbation strength range, a 
set of twenty numerical prediction realizations were performed with V a set first to 0.125 
m/s and then twenty more realizations with V a set to 0.25 m/s. These strength parameters 
produced maximum sound speed perturbations of approximately 6 m/s and 12 m/s respec- 
tively. 

For each of the twenty range-dependent realizations with V 0 set to 0.125 m/s, the 
phase a.j was randomly assigned. However to allow direct comparison between data sets, 
the same set of random phases was used for the twenty range-dependent realizations with 
V Q set to 0.25 m/s. Table 2 lists the assigned phase values used for each of the four 
baroclinic modes in each realization. Additionally, the location of the local sound speed 
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axis at the initial range step was directly influenced by the phase of the random mesoscale 
structure and the perturbation strength. Table 2 also provides the resulting depth of the 
local sound speed axis and the corresponding minimum sound speed at this local axis for 
each perturbation strength. In all realizations, the assigned mesoscale sound speed 
structure is assumed “frozen” in time over the pulse propagation period. Small-scale 
fluctuations such as those caused by internal waves are not considered. 

A typical realization of the sound-speed structure for V o =0.125 m/s is shown in 
Fig. 6. The perturbations in this particular realization range from approximately -5.3 to 
6. 1 m/s relative to the background sound speed profile. Figure 6 also provides the sound- 
speed perturbation envelope referenced to the range- independent canonical profile. Figure 
7 provides a similar illustration of the sound-speed structure for the case V o =0.25 m/s. In 
this case, the perturbations range from -10.7 to 12.2 m/s. For both perturbation strengths, 
the perturbation is essentially confined to the upper half of the water column. Note that 
the range-dependence of these perturbations is sinusoidal. Although the sound speed 
perturbation will average to zero over long range, the travel time perturbation is non-linear 
and will not average to zero. 

B. ACOUSTIC PROPAGATION MODELING 

Numerical predictions of the acoustic propagation were obtained from a modified 
version of the UMPE model (Smith and Tappert, 1993). UMPE is a broadband, range- 
dependent parabolic equation model which implements the split-step Fourier (SSF) 
method described previously in Chapt. II. It was modified to include a subroutine for 
generating the random mesoscale perturbation and superimposing these perturbations onto 
the background sound speed profile. Input run parameters used for acoustic propagation 
modeling are summarized by Table 3. 

A continuous wave pulse, wide-angle acoustic source function was used to 
generate the starting field for multiple frequencies. The solutions are then inverse Fourier 
transformed to yield travel times. The modeled source has a 100 Hz center frequency, 
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Randomly assigned phase value for 
baroclinic mode number 


V 0 =O.125 


V 0 =0.25 


Run 


1 


2 


3 


4 


Axis 

Depth 

(m) 


Sound 

Speed 

(m/s) 


Axis 

Depth 

(m) 


Sound 

Speed 

(m/s) 


1 


0.4437 


2.0744 


0.0500 


3.5025 


1039 


1491 


1100 


1492 


2 


2.9363 


6.1271 


0.9400 


1.2410 


836 


1489 


695 


1486 


3 


2.4603 


2.2026 


6.1550 


1.2024 


898 


1488 


845 


1486 


4 


1.4359 


3.2068 


3.1620 


1.8819 


1074 


1490 


1118 


1490 


5 


5.6134 


1.9929 


5.1243 


3.3137 


1065 


1491 


1136 


1492 


6 


5.6213 


6.0501 


5.7924 


3.3372 


1012 


1492 


1083 


1494 


7 


1.0054 


0.1503 


5.9521 


1.7476 


962 


1491 


872 


1492 


8 


2.6148 


2.1163 


2.1666 


1.5559 


971 


1488 


951 


1487 


9 


2.6314 


1.6137 


0.2810 


3.0539 


1050 


1491 


704 


1486 


10 


1.0813 


5.4430 


4.9802 


6.2065 


1033 


1490 


986 


1491 


11 


2.5339 


2.9236 


1.3976 


0.6285 


945 


1488 


916 


1486 


12 


0.8081 


3.6260 


3.2656 


0.7081 


945 


1488 


1136 


1490 


13 


1.5948 


2.9414 


1.8233 


3.9961 


1024 


1489 


1065 


1489 


14 


0.6205 


6.0390 


0.7910 


5.8630 


1042 


1492 


986 


1493 


15 


0.4336 


0.7553 


1.8742 


4.2228 


1042 


1491 


1197 


1493 


16 


5.3948 


5.1485 


5.0633 


4.8832 


1138 


1491 


1047 


1492 


17 


1.4660 


0.9852 


6.0588 


6.1911 


1006 


1490 


880 


1490 


18 


0.2466 


1.7549 


2.9578 


2.5819 


989 


1488 


124T 


1492 


19 


3.2380 


3.7665 


0.2168 


0.3478 


998 


1490 


854 


1484 


20 


1.0107 


3.6465 


2.8894 


6.2395 


1050 


1492 


1109 


1490 



Table 2. Assigned baroclinic mode phase values with resulting local sound speed axis 
depths and axis sound speeds at the initial range step. 
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Sound speed perturbation envelope, Vo = 0.125 m/s 




Figure 6. Sound speed perturbation field realization for V0=0.125 m/s. The upper panel 
displays the sound speed perturbation field relative to the unperturbed background sound 
speed field. The lower panel compares the background sound speed profile with the 
envelope of the perturbed profiles. 
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Typical perturbation sound-speed structure, Vo=,0.25 




Sound speed perturbation envelope, Vo = 0.25 m/s 




Sound speed (m/s) 



Figure 7. Sound speed perturbation field realization for V 0 =0.25 m/s. The upper panel 
displays the sound speed perturbation field relative to the unperturbed background sound 
speed field. The lower panel compares the background sound speed profile with the 
envelope of the perturbed profiles. 
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Parameter 


Value 


Type of PE approximation 


Wide-angle 


Reference sound speed, c 0 


1490.0 m/s 


Number of points in depth 


2048 


Computational depth 


9000.0 m 


Depth mesh size 


8.8 m 


Range step size 


50.0 m 


Maximum range 


1000.0 km 


Center frequency 


100.0 Hz 


Number of frequencies 


512 


Frequency bandwidth 


100.0 Hz 


Frequency bin spacing 


0.195 Hz 


Source depth 


Varied, placed at local sound speed axis, 
which varied up to ±200 m from 
background sound speed axis at 1000 m. 


Source level 


200 dB re 1 jx Pa 


Source spectrum 


Hanning amplitude window 


Angle of the beam center 


0.00 deg 


RMS half-width angle 


35.00 deg 


Taper 


15.00 deg 



Table 3. Acoustic propagation model (UMPE) input run parameters for deep ocean 
environment with mesoscale structure. 
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100 Hz bandwidth sampled at 512 discrete frequencies, and 200 dB re lflPa source level. 
The UMPE code generates the frequency bin spacing and assumes a Hanning amplitude 
window for the source spectrum. The pulse length is not an input parameter. A wide- 
angle source starting function (Thomson and Chapman, 1983) was used with a D/E angle 

of 0° and a half-width of 35°. A smooth filter was applied to the outer 15° to reduce the 
influence of sidelobes. To ensure the lowest modes were excited in each realization, the 
source was placed at the local sound speed axis, the depth of which varied by as much as 
200 meters from the assumed background axis depending on the specific perturbation 
field. 

To numerically predict the propagation field, a computational depth of 9000 m 
was used with a total transform size of 2048, thus providing for a depth mesh of 8.8 m. 
The projected field was computed out to a range of 1000 km using a range step of 50 m. 
Twenty range-dependent realizations for each perturbation strength V 0 were computed 
using the randomly-phased baroclinic modes as per Table 2. Next, a range -averaged 
sound speed profile was computed for each realization using the range-dependent sound 
speed profiles at 1 km steps up to the range of interest. To investigate the bias in travel 
time resulting from the mesoscale structure, forty additional realizations for each pertur- 
bation strength were then computed. Twenty each realizations were computed out to 500 
km and 1000 km ranges using the respective range-averaged sound speed profile. Note, 
that while this use of the range-averaged sound speed profile will result in a range- 
independent environment, the range-averaged sound speed profile may differ slightly from 
the background Munk canonical profile (no perturbation). From the computed fields, the 
solutions are inverse Fourier transformed to yield travel times. 

A typical plot of the computed broadband transmission loss field for the range- 
dependent case with v 0=0 .125 m/s is provided in Fig. 8. In computing the FFTs, the time 

domain is heterodyned around the value t Q = /?/c 0 . The time scale axis then uses 
“reduced time” or (t-t 0 ) for the arrival time. Some wrap-around of the signal is observed 
due to the finite transform size. Note that the arrival structure even for the latest arrival is 
well-defined, suggesting an adiabatic nature of the propagation. 



51 



Depth (m) Depth (m) 



60 



120 



. 

Transmission Loss (dB re 1m) 



0 

1000 

2000 

3000 

4000 

5000 




- 2.5 - 2.0 - 1.5 - 1.0 - 0.5 0.0 0.5 1.0 1.5 2.0 2.5 



Reduced Time (sec) 



0 



1000 



2000 




3000 

- 0.5 



- 0.4 



- 0.3 



- 0.2 



- 0.1 0.0 0.1 
Reduced Time (sec) 



Figure 8. Transmission loss arrival time structure at range 1000 km for V 0 =0.125 m/s. 
Upper panel shows total structure of the range-dependent realization. Lower panel shows 
expanded view of axial arrivals. 
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To further assess the relative adiabatic nature of the modeled environment, a modal 
decomposition was performed on a single frequency (100 Hz) range-dependent realization 
(#7) using the second order approximate wide angle mode functions described in Chapt. 
III. Results are provided by Fig. 9. At perturbation strength V^pO.125 m/s, there is only a 
limited amount of mode coupling. However, as the perturbation strength is increased to 
Vo=0 .25 m/s, mode coupling starts to increase. 

Representative arrivals at the 1000 m depth and 500 km range are depicted in 
Fig. 10. In the upper panel, the arrival structure of a range-dependent realization for 
V fy=0.25 m/s is compared with its range- averaged realization. The lower panel shows an 
expanded version of the first and last arrivals. Figure 1 1 provides similar illustrations for 
the 1000 km range. In all but one of the range-dependent realizations computed for both 
Vq= 0.125 m/s and V(j= 0.25 m/s, the latest peak arrival at the 1000 km range was delayed 
from that of the range-averaged case. 

C. DATA ANALYSIS AND DISCUSSION 

For this modeled environment, we are interested in characterizing the overall 
spreading in the arrival structure envelope of the range-dependent realization relative to its 
range-averaged counterpart. In particular, we wish to determine if there is a statistically 
significant bias in the near-axial arrival structure. The data analysis will refer to the data 
in four sets based on the combinations of range (500 km, 1000 km) and perturbation 
strength (0.125 m/s, 0.25 m/s). 

The first step in analyzing the data was to measure the arrival envelope lengths. 
Employing the technique of using arrival times as the time instants corresponding to the 
significant maxima (peaks) of the arrival pattern (Athanssoulis and Skarsoulis, 1995), the 
arrival envelope length was determined by measuring the time difference between peaks 
of the first and last arrivals. Only peaks with less than 100 dB re 1 m transmission loss 
were considered part of the arrival structure. In those cases where the first arrival 
“wrapped around” due to the finite transform window, the relative time was adjusted 
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Vo=0 




Vo=0.25 




Figure 9. Representative modal decomposition of computed single frequency (100 Hz) 
pressure fields. Amplitudes for the ten modes having the largest amplitudes are plotted. 
Only the three modes with the largest amplitudes are identified. 

54 



TL, Mesoscale Vo=0.25, Depth=1000m, Range=500 km 
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Figure 10. Transmission loss arrival time structures at range 500 km and depth 1000 m. 
Upper panel shows total structure for range-independent and range-dependent realization 
with V0=O.25 m/s. Lower panels show expanded view of the first and last arrivals. 
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TL, Mesoscale Vo=0.25, Depth=1000m, Range=1000 km 




Expanded view of first arrival Expanded view of last arrival 





Figure 11. Transmission loss arrival time structures at range 1000 km and depth 1000 m. 
Upper panel shows total structure for range-independent and range-dependent realization 
with V0=O.25 m/s. Lower panels show expanded view of the first and last arrivals. 
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accordingly. Histograms showing the respective number of realizations in each 10 msec 
time bin are provided by Figs. 12 and 13 for the 500 km and 1000 km ranges, respectively. 
Note that as the perturbation strength is increased, there is a significant spreading in the 
arrival envelope lengths. Additionally, the mean range-dependent arrival envelope length 
is longer than its corresponding mean range-averaged length for all four sets of data 
plotted. 

The arrival envelope spread is computed by taking the difference between the 
range-dependent and associated range-averaged arrival envelope lengths. A positive 
spread indicates that the range-dependent envelope length is greater than its associated 
range-averaged counterpart. Again, histograms with 10 msec bins were constructed and 
are provided in Figs. 14 and 15 for the 500 km and 1000 km ranges, respectively. Note 
that while the mean for all four data is positive, there are several individual realizations 
with negative spreads. However, the arrival envelope spread becomes increasingly 
positive with increasing range and perturbation strength. 

Next, the time difference between the range-dependent and associated range- 
averaged first peak arrivals was measured. This same time difference was also measured 
for the last peak arrivals. In this case, we refer to the time difference as the bias, where 

X BIAS = X RD~ X ra‘ - ( 4 - 6 ) 

In Eq. (4.6), the subscript RD denotes range-dependent and RA denotes range-averaged. 
Figures 16 and 17 provide histograms of the bias in the first and last peak arrivals at the 
500 km and 1000 km ranges, respectively. 

The mean bias for the near-axial (last) arrivals is generally positive, and increases 
with increasing perturbation strength and range. In particular, the mean bias for the last 
peak arrival increases by over 100 msec when the perturbation strength is increased from 
Vq= 0.125 m/s to .25 m/s. A possible physical explanation for this cold bias is that the 
acoustic path is continually adjusting toward the lowest sound speed, thus resulting in an 
actual sound speed along the path which is lower than the comparable range-average. The 
mean bias for the higher modes is significandy less indicating that these modes are more 
stable in terms of their arrival times. In some realizations, the earlier arrivals, those 
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Figure 12. Arrival envelope lengths at 500 km range. 
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Figure 13. Arrival envelope lengths at 1000 km range. 
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Figure 14. Arrival envelope spread at 500 km range. 
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Figure 15. Arrival envelope spread at 1000 km range. 
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Figure 16. Travel time bias in first and last arrival peaks at 500 km range. 
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Travel time bias in first and last arrival peaks at 1000 km range. 
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corresponding to higher acoustic modes, actually arrive slightly earlier when the 
perturbations are present. The data may then suggest that the total time envelope of 
arrivals (i.e. time between earliest and last arrivals) always tends to increase. In all cases, 
however, the influence is most pronounced in the delayed arrival of the lowest modes. 

In closing out this chapter, there are four main points which are suggested by the 

data. 

1) Mesoscale features can significantly perturb the late arrivals (corresponding to the 
lowest order acoustic modes) and produce an overall spreading in the arrival structure 
even when propagation is adiabatic. 

2) The bias in the arrival times of the lowest modes increases with increasing perturba- 
tion strength. 

3) The size of bias variability is not directly related to range or perturbation strength. 

4) Steeper, earlier arrivals are relatively stable in the presence of mesoscale structure. 
The arrival data statistics are summarized by Table 4. 
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Range 


Range=500 km 


Ran ge= 1000 km 


Strength 




Mean 


Std Dev 


Mean 


Std Dev 


Vo= 

0.125 m/s 


First peak 
bias 


-3.1 ms 


8.4 ms 


16.8 ms 


19.6 ms 


Last peak 
bias 


25.5 ms 


35.4 ms 


88.4 ms 


74.1 ms 




Envelope 

Spread 


28.6 ms 


34.6 ms 


71.6 ms 


56.5 ms 


Vo= 

0.25 m/s 


First peak 
bias 


-10.0 ms 


23.8 ms 


9.4 ms 


8.9 ms 


Last peak 
bias 


89.6 ms 


79.2 ms 


205.0 ms 


55.8 ms 




Envelope 

Spread 


99.6 ms 


77.9 ms 


195.8 ms 


54.7 ms 



Table 4. Arrival statistics for deep ocean with random mesoscale structure. 
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V. TRAVEL TIME VARIABILITY DUE TO THE CALIFORNIA CURRENT 

FRONT 



In this chapter, the effects of a fluctuating frontal region created by the California 
Current on modal acoustic travel time variability are investigated using a state-of-the-art 
global ocean circulation model to provide input environmental data. The California 
Current is the major eastern boundary current of the North Pacific. It is 500-1000 km 
wide and has been found to be highly variable in both space and time (Mooers and 
Robinson, 1984). The instantaneous current consists of intense meandering current 
filaments intermingled with synoptic-mesoscale eddies. In general, cool water upwelled 
near shore is entrained and then subsequently advected 100 km or more offshore by the 
variously oriented jetlike strong current filaments and eddies. Additionally, the jets and 
eddies can also advect warm offshore waters onshore to the coastal ocean. Through 
satellite imagery and oceanographic studies, it has been shown that this jet and eddy 
system can change substantially on a weekly time scale. 

Recently there has been work along a similar theme (Dushaw, 1997) which was 
motivated by the ATOC project. Dushaw’s concern was that, as described previously, 
variability in the acoustic pulse travel time acts as geophysical “noise” in problems such 
as accurately assessing the large-scale North Pacific heat budget. An accurate assessment 
may be extremely difficult for acoustic paths passing through this region if the effect of 
variability in the California Current is severe enough. Dushaw examined the effects of the 
California Current on synthetic acoustic transmissions made from the ATOC Pioneer 
Seamount acoustic source (shown on Fig. 18) to a fictional receiver 500 km away along a 
path which is nearly perpendicular to the California coast. He calculated the variation of 
ray travel times using environmental input data interpolated from the California Cooper- 
ative Fisheries Investigations (CalCOFI) time series data. However, there are shortfalls in 
using the CalCOFI data which are pointed out by Dushaw. In particular, the CalCOFI 
data sampling is irregular in time, depth, and horizontal sampling. At best, temporal 
sampling is provided at some of the CalCOFI grid stations on a monthly basis, but the 
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entire grid is rarely surveyed during any particular month or year. The sporadic monthly 
sampling does not capture the complete variability known to exist in this region from 
satellite imagery. Furthermore, vertical sampling is limited to the upper 500 to 1000 m of 
the water column, dependent on the observation year. Accordingly, Dushaw relied on a 
limited number of sound speed profiles and extensive interpolation to generate his range- 
dependent sound speed environment, and as discussed in his results, his study could not 
address any sound speed variability which may occur below 1000 m. Given these limita- 
tions, and after accounting for the assumed general North Pacific variability, his 
conclusion was that a conservative upper bound on ray travel time variability specifically 
due to California Current variability is ±50 ms for the 1000-5000 km ATOC transmis- 
sions. 

In an effort to capture more completely the variability of the California Current 
System and acoustic transmissions through a small region of it, this investigation uses the 
global Parallel Ocean Circulation Model (POCM) (Semtner and Chervin, 1992), also 
known as the Semtner and Chervin model, to obtain realistic environmental input data 
spanning a one-year time period over a 474 km acoustic path. In that this was a first-time 
use of POCM with the UMPE acoustic propagation model, it was necessary to develop 
computer code to interface the two models. Then, using the coupled ocean-acoustic 
model, the temporal, spatial and seasonal variability in the individual modal arrivals of the 
first thirty modes is assessed. The mesoscale bias variability is also examined by 
comparing the full-field peak arrival times for the range-averaged environment to that of 
the range-dependent environment at the 400 km and receiver ranges. 

The research here is presented in a number of sections. First the ocean model used 
for environmental data input is described. Next the method of reformatting the input 
environmental data into a usable format for UMPE and description of a typical sound 
speed field are provided. This is followed by a discussion of the acoustic modeling 
method and modal variability analysis. Finally, a discussion of results with general 
conclusions is presented. 



68 



A. 



OCEAN MODEL 



Input environmental data is obtained from the global Parallel Ocean Circulation 
Model (POCM), also known as the Semtner and Chervin model. Semtner and Chervin 
(1992) describe the basic POCM model formulation, and Stammer et. al. (1996) describe 
subsequent improvements to the POCM model. It was this improved model that was used 
for obtaining input environmental data for the coupled ocean-acoustic modeling described 
herein. Key features of the improved model are described below. 

The POCM model is capable of resolving mesoscale ocean features and strong 
currents both spatially and temporally. Initially, a 1/2-degree model initialized with 
Levitus temperature and salinity data was spun up for 33 years. Instantaneous fields were 
then interpolated to a Mercator grid size of 0.4 degree in longitude, thus producing square 
grids between the equator and 75 degrees latitude. This provides a nominal lateral 
resolution of 1/4 degree (-25 km). The simulation was then resumed, starting with the 
recorded winds of January 1985 providing the forcing function. Up to 20 vertical layers 
are provided at each grid point. In that the greatest variability is observed near the ocean 
surface, the spacing between vertical layers is approximately 25 m near the surface and 
then the layer spacing gradually increases with ocean depth. Coastlines and bathymetry 
are prescribed at the local model resolution, and are not smoothed further. 

Model data is available for numerous years, at a time interval of three days. With 
realistic external surface boundary forcing, including daily wind-stress fields and sea 
surface heat fluxes, the POCM simulated ocean environment successfully models the 
large scale variability as compared with oceanographic observations. However, the 
simulated amplitudes of variability are reported to be low by about a factor of 2 to 4 over a 
broad spectral range. This discrepancy in energy is greatest at high frequencies and 
wavenumbers, diminishing to a factor of two at the lowest frequencies and wavenumbers 
(Stammer et al., 1996). 
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B. ACOUSTIC PATH DESCRIPTION AND INPUT DATA PREPARATION 



For this particular investigation, the chosen acoustic path is between a moored 
SOFAR source designated SS3 and a SOSUS hydrophone array at Point Sur, California 
for an approximate path length of 474.6 km. The source and receiver locations are shown 
in Fig. 18. The Naval Postgraduate School (NPS) sound sources designated SSI, SS2, 
SS3 in this figure were specifically deployed to study the California Current System by 
tracking oceanographic floats (Garfield et. al., 1997). The source provides a linear 
frequency modulated chirp with a center frequency of 260.1 Hz and a swept bandwidth of 
1.523 Hz, has a source level of 183 dB re 1 |iPa, and is moored at 642 m depth. Source 
signal monitoring is done with the Pt. Sur SOSUS hydrophone array (C.-S. Chiu, pers. 




Figure 18. West coast sound sources. NPS sound sources designated SS 1 , SS2, SS3 were 
specifically deployed to study the California Current System by tracking oceanographic 
floats. SOSUS array is located near Pt. Sur. 
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comm., 1997). The Point Sur array is at the ocean bottom, a depth of 1359 meters. Actual 
data is available from this source/receiver combination for future comparison with 
numerical predictions obtained using POCM environmental input data. Available source 
and receiver characteristics are summarized in Tables 5 and 6, respectively. 

To assess the acoustic stability of arrivals for this source/receiver combination, the 
first step was to obtain data sets of temperature and salinity spanning a one-year time 
period (122 data sets, model year 1995) from the POCM model. Bathymetry data was 
obtained from the National Oceanographic Data Center (NODC) database DBDB5 which 



Parameter 


Value 


Mooring location 


37° 06.55’N 127° 34.60’W 


Source depth 


642 m 


Water depth 


4752 m 


Center frequency 


260.1 Hz 


Signal 


LFM Chirp 


Swept Bandwidth 


1.523 Hz 


Pulse length 


0.6 sec 


Source Level 


183 dB re 1 |iPa 



Table 5. Naval Postgraduate School SS3 source characteristics 



Parameter 


Value 


Location 


36° 17.87’N 122° 22.58’W 


Type 


Single hydrophone 


Receiver depth 


1359 m 


Bottom depth 


1359 m 


Distance from SS3 source 


476.1 km 



Table 6. Pt. Sur receiver characteristics 
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has a resolution of 0° 05’. The next step was to interpolate the temperature, salinity, and 
bathymetry data from the ocean model/database grid points onto grid points along a 2-D 
great circle path between the chosen source-receiver combination. Horizontal 
interpolations were done linearly, whereas interpolations between vertical depth layers 
were done using a cubic spline. The interpolated temperature and salinity data was then 
converted into sound speed using the empirical formula (Mackenzie, 1981) 

c(z,S,t)= 1 449.2 + 46/- 5.5^ + 0.29^ + (1.34-0.1r)(5-35)+ 0.0 16z (5.1) 

where c is sound speed in m/s, 5 is the salinity in parts per thousand (ppt), z the depth in 
meters, and /=7710 where Tis in degrees Celsius. 

A typical sound speed profile interpolated from the POCM output consists of a 
surface duct with its local sound speed minimum around 200 m depth and a main water 
column channel with an axial depth around 1000 m. To illustrate the annual sound speed 
variability along the chosen acoustic path, the yearly mean range-averaged sound speed 
profile and its envelope are provided in Fig. 19. This figure also plots the full ensemble of 
range-averaged sound-speed profiles. Note that the POCM output suggests that the 
surface duct persists throughout the annual cycle. This is not considered realistic (C.-S. 
Chiu, pers. comm., 1997). Subsequently, this chapter will focus on modal travel time 
variability for modes which remain in the main water column channel. 

The sound speed perturbation fields for POCM model dates 05 Jan 95 and 01 Jul 
95 are depicted by Figs. 20 and 21, respectively, where the sound speed perturbation is 
defined as 

Sc(r, z) = c(r,z) ~c RA (z) (5.2) 

and c RA (z) is the reference range-averaged sound speed profile specific to that data set. 

For the January perturbation field, the sound speed perturbations range from about ±3.5 
m/s and are essentially limited to the upper 2000 meters in the water column. The July 
perturbation field shows much stronger gradients, with amplitudes ranging from about 
±6.0 m/s and extending deeper to near 2500 meters depth. Throughout the year, the 
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Figure 19. Annual range-averaged sound speed profiles as interpolated from POCM 
output for modeled environment. Upper panel shows the mean annual range-averaged 
sound speed profile and envelope of range-averaged sound speed profiles. Lower panel 
illustrates the variation in range-averaged sound speed profiles over the annual cycle. 
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Figure 20. Sound speed perturbation field for POCM model date 05 Jan 95. Top left panel 
shows reference range-averaged sound speed profile. Top right panel shows 1 1 equally- 
spaced sound-speed profiles as a function of range. Abscissa scale applies only to first 



profile. Middle panel shows perturbation field for entire water column and bathymetry 
along track. Lower panel provides enlarged view of field in the upper 400 m. 
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Figure 21. Sound speed perturbation field for POCM model date 01 Jul 95. Top left panel 
shows reference range-averaged sound speed profile. Top right panel shows 1 1 equally- 
spaced sound-speed profiles as a function of range. Abscissa scale applies only to first 
profile. Middle panel shows perturbation field for entire water column and bathymetry 
along track. Lower panel provides enlarged view of field in the upper 400 m. 
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surface duct is seen to be much more pronounced in the region near the receiver relative to 
that near the sound source. 

C. ACOUSTIC MODELING 

The range-dependent sound speed environments computed using the 122 POCM 
datasets are individually input into the UMPE acoustic propagation model to compute the 
full acoustic pressure field and the arrival structures at specific ranges of interest. In 
addition, UMPE runs using a range-averaged sound speed profile as environmental input 
were completed for 25 data sets (every 15 model days) to support analysis of the 
mesoscale bias variability. Range- averaged sound speed profiles were computed for each 
data set by averaging the range-dependent sound speed profiles (water column sound 
speeds only) at 1 km steps up to the range of interest. Input parameters for the UMPE 
acoustic propagation modeling are summarized in Table 7. Typical range-dependent 
arrival structures across the full water column at the receiver range and two intermediary 
ranges are presented in Fig. 22. Figure 23 displays a typical full-field arrival structure in 
terms of transmission loss specifically at the receiver depth and range. The variability in 
these full-field arrival structures at the receiver depth and range are demonstrated by Figs. 
24 and 25. 

To determine the degree of mode coupling as a function of environmental 
parameters and bathymetry, pressure field modal decompositions were completed for 
range-averaged and range-dependent environments, with and without bathymetry. The 
standard PE approximation was used to compute the pressure field used for these modal 
decompositions, so that KRAKEN standard normal modes might be employed as the 
proper basis set. Typical results are plotted in Fig. 26. For the modes analyzed, this figure 
clearly illustrates that bathymetry plays no role in mode coupling until approximately the 
last 20 km where the ocean bottom shows a relatively steep rise to depths of less than 3000 
m as it nears the Pt. Sur coastline. Environmental parameters which influence sound 
speed appear to have the greatest impact on mode coupling in the 200-300 km range 
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Parameter 


Value 


Type of PE approximation 


Wide-angle 


Reference sound speed, c 0 


1495.2 m/s 


Number of points in depth 


8192 


Computational depth 


6500.0 m 


Depth mesh size 


1.6 m 


Range step size 


5.0 m 


Maximum range 


474.61 km 


Center frequency 


260.1 Hz 


Number of frequencies 


64 


Frequency bandwidth 


1.6 Hz 


Frequency bin spacing 


0.025 Hz 


Source depth 


642.0 m 


Source level 


183 dBre 1 ^iPa 


Source spectrum 


Hanning amplitude window 


Angle of the beam center 


0.00 deg 


RMS half-width angle 


35.00 deg 


Taper 


10.00 deg 



Table 7. Acoustic propagation model (UMPE) input run parameters for west coast 
environment. 

segment, and again between 450 and 474 km range. In this latter range segment, surface 
ducting and bathymetry effects become important to the specific modes being analyzed. 
As will be shown in the next section, some of the local normal modes shift from a “deep 
water” mode at ranges of 200 and 400 km to surface duct modes at the receiver range. 
Due to concerns regarding the validity of surface layer sound speeds as interpolated from 
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Figure 22. Typical range-dependent arrival structures at 200, 400, and 474 km ranges. 
Reduced time is referenced to 133.761, 267.523, and 317.081 seconds, respectively. 
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Range=474.6km, depth= 1359m 




Figure 23. Plot of typical arrival structure at receiver depth (model date 31Jul95). 




Figure 24. Full pressure field arrival structure at receiver over one model year period. 

325 r 




Figure 25. Centroid and energy distribution of arrival structure at receiver over one year. 
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RANGE-AVG SSP w/BATHYMETRY 






Figure 26. Modal coefficient amplitudes for a single frequency run. Upper panel is for the 
range-averaged sound speed profile with bathymetry. Middle panel assumes flat bottom 
with range-dependent sound speed profile. Lower panel incorporates both the range 
dependent sound speed profile and bottom bathymetry. 
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POCM output, modal travel time variability was assessed only at the 400 km range where 
the first thirty modes remained in the main water column channel. 

D. MODAL ANALYSIS 

In this section the spatial and temporal variability of the first thirty local modes is 
assessed. In addition, the mesoscale bias variability is examined by comparing the full- 
field peak arrival times for the range-dependent environment to that of the corresponding 
range-averaged environment. To examine the modal variability, the WAPE pressure field 
at selected ranges was decomposed into local normal modes using the numerical code 
KRAKEN (Porter, 1995). As discussed in Chapt. HI, there is a slight, but acceptable, 
error introduced to the modal travel time by using standard normal modes for decom- 
posing the WAPE field. Using inverse Fourier techniques described in Chapt. II, the 
individual modal arrival times were computed at the 400 km range. Note that due to mode 
coupling, there is a complex relationship between modal travel time variability and sound 
speed perturbation. 

The spatial variability in selected local mode functions at the 200 km, 400 km, and 
receiver ranges are depicted by waterfall displays in Figs. 27, 28, and 29, respectively. 

At the 200 and 400 km ranges, modes 1-30 reside entirely in the main water column 
channel, with upper and lower turning points confined between 400 m and 1950 m depths. 
At the 200 km range, the mode shapes show very little variability over the annual cycle. 
At the 400 km range, there is a slight but detectable upward shift in the upper and lower 
turning points during the summer season. The greatest spatial variability in the mode 
functions occurs at the receiver range. At this range, modes 23-30 alternately appear as 
surface duct (modes 1 and 2) and water column channel modes dependent on the specific 
sound speed environment and frequency. For simplicity in further discussion, these 
modes will be referred to simply as “surface duct” modes. The axis of the surface duct 
clearly shifts over time, as is evident in the mode 30 plot. Again, the modal upper turning 
points generally tend to shift upwards during the summer season. In addition, it is noted 
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Figure 27. Annual spatial variability of selected mode functions at 200 km range. 
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Figure 28. Annual spatial variability of selected mode functions at 400 km range. 
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Figure 29. Annual spatial variability of selected mode functions at receiver range. For 
modes 23 and above, mode numbers refer to deep channel mode number (rather than 
surface duct mode number). 
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that modes 8 and higher are bottom interacting at the receiver range. 

To assess the temporal variability for each of the first 30 modes, the following 
statistics were computed at the 400 km range: 

1) ensemble (annual) mean for the energy centroid arrival time, and the 
corresponding ensemble mean energy distribution (arrival width), 

2) variance in the ensemble mean centroid arrival time, and 

3) variance in the ensemble mean energy distribution, 

where the ensemble consisted of the 122 range-dependent data sets. Results are presented 
graphically for selected modes in Figs. 30 and 31, and numerically for the first thirty 
modes in Table 8. These statistics show that arrival time (centroid) variability is limited to 
±52 msecs. Modes 12 through 18 are the most stable, with their variability limited to 
±17 msecs. The mean energy distribution (spread in modal peak arrival) also shows little 
variability, in general less than ±6 msecs, over the annual cycle for these lower water 
column modes with the exception of modes 1, 2, 3, and 30 which vary up to ±276 msec. 
However, it is also noted that for many of the computations, modal amplitudes for modes 
1 and 2 were so low that sidelobes in the FFT process might be affecting the results. 

To assess the seasonal variability, the frequency spectra of the modal arrival times 
was computed. Figure 32 plots the modal arrival time frequency spectra for selected 
modes. The dominant frequency of variability for all modes examined is seasonal (1 
cycle/year). Higher harmonics also appear for modes 1, 2, 29, and 30. 

The mesoscale bias variability is examined by comparing the full-field peak arrival 
time for the range-dependent environment to that of the corresponding range-averaged 
environment at the 400 km and receiver ranges. In determining the range-averaged 
environment, only sound speeds in the water column were averaged. With bias defined as 
the difference between the range-dependent and the range-averaged travel times, 
x bias = x rd ~ x ra * computations included determining 

1) the mean bias over 25 samples, and 

2) the associated standard deviation for the mean bias. 



85 



MODE 1 MUOt a 




MODE 13 MODE 17 




MODE 23 MODE 30 




Figure 30. Annual variability in centroid of select modal arrivals at the 400 km range. 
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Figure 3 1 . Temporal variability in first 30 modes at the 400 km range. Upper panel left 
plots the modal ensemble mean centroid arrival time with mean energy distribution. 
Upper right panel shows variance in the modal ensemble mean arrival time. Lower left 
panel shows the variance in the modal ensemble mean energy distribution. Lower right 
panel combines information presented in first three to show the ensemble variance in the 
modal arrivals and energy distribution. 



87 



Mode # 


Ensemble mean 
centroid arrival 
(secs) 


Std deviation of 
ensemble mean 
centroid arrival (ms) 


Ensemble mean 
energy 

distribution (ms) 


Std deviation of 
ensemble mean energy 
distribution (ms) 


1 


269.98 


51.7 


528.2 


275.3 


2 


269.99 


36.4 


446.4 


95.7 


3 


269.99 


27.0 


399.6 


25.6 


4 


269.99 


20.7 


389.6 


6.0 


5 


269.99 


19.8 


388.1 


1.5 


6 


269.98 


18.6 


388.2 


0.4 


7 


269.98 


21.2 


388.2 


0.7 


8 


269.97 


17.7 


388.2 


0.4 


9 


269.96 


17.5 


388.4 


0.8 


10 


269.95 


18.9 


387.5 


0.5 


11 


269.95 


18.4 


387.7 


0.3 


12 


269.94 


16.7 


387.7 


0.1 


13 


269.93 


16.8 


387.5 


0.1 


14 


269.93 


16.4 


387.5 


0.1 


15 


269.92 


16.0 


387.2 


0.1 


16 


269.92 


15.7 


387.2 


0.1 


17 


269.92 


15.9 


387.6 


0.1 


18 


269.92 


17.0 


387.6 


0.1 


19 


269.93 


19.2 


387.8 


0.2 


20 


269.93 


21.9 


387.7 


2.1 


21 


269.94 


21.2 


387.4 


0.7 


22 


269.94 


21.7 


388.0 


1.5 


23 


269.94 


18.9 


387.5 


0.2 


24 


269.93 


18.6 


387.9 


0.1 


25 


269.93 


18.6 


388.0 


0.2 


26 


269.93 


19.7 


387.8 


0.1 


27 


269.92 


23.0 


388.3 


0.5 


28 


269.92 


22.2 


388.2 


0.6 


29 


269.91 


30.0 


388.5 


3.4 


30 


269.90 


27.8 


405.3 


188.6 



Table 8. Ensemble statistics for modal arrivals at the 400 km range. 
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Figure 32. Frequency spectra of modal arrival times at the 400 km range. 
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For these calculations, the full-field arrival structure is dominated by a single peak arrival, 
as is evident in Fig. 23. This is partially due to the small bandwidth of the source 
spectrum. In addition, the source depth was generally several hundred meters above the 
deep sound channel axis. Thus, it is the higher modes (steeper propagation angles) which 
are predominantly excited. For the bias analysis, only this single peak was analyzed. 

Histograms of the sampled bias at the two ranges are provided by Fig. 33. Bin 
spacing is 156 ms. In contrast to the generally cold bias results shown for the open ocean 
environment analyzed in Chapt. TV, the mean bias for the 400 km range appears to have a 
Gaussian-like distribution about zero bias. As previously mentioned, for this environment 
and source placement, only higher modes are predominantly excited. As suggested by the 
open ocean environment analysis, steeper arrivals tend to be more stable in the presence of 
mesoscale structure. Thus, this bias result does not contradict the general findings of 
Chapt. IV. At the 474 km range, the distribution indicates a tendency towards a negative 
(warm) bias. One possible factor contributing to this result for this acoustic path is the 
effect of bathymetry. With the steep rise in bathymetry near the receiver, energy is 
stripped out of the highest (and generally faster) modes and subsequently transferred into 
the ocean bottom or coupled into lower modes. Due to the range-averaging of the sound 
speed structure, the range-averaged ocean environment near the receiver is warmer than 
the range-dependent case. Any coupling would then shift energy into faster modes than 
predicted in the range-dependent ocean, thereby causing a general warm bias trend. Thus, 
this environment would suggest that mesoscale travel time bias through a region with 
significant bathymetry influences is not a simple function of range. Finally, the large bias 
variability, as indicated by a standard deviation in the mean annual bias in excess of 500 
ms at each of the two ranges, suggests that the ability to accurately interpret travel time 
trends by assuming a range-averaged sound speed profile may be limited. 
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Peak arrival bias at 400 km range 




Peak arrival bias at 474 km range 




Figure 33. Bias in full-field peak arrival. Upper panel shows bias at the 400 km range, 
while lower panel shows bias at 474 km receiver range. Statistics based on a sample size 
of 25 with temporal sampling every 15 days. 
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E. 



SUMMARY 



In this chapter, acoustic modal travel time variability through a region of the 
California Current was investigated using a coupled ocean-acoustic model. Basic results 
from this investigation are interpreted with the caveat that the simulated amplitudes of 
mesoscale variability in the oceanographic model are known to be low by about a factor of 
2 to 4 over a broad spectral range. With that in mind, we have seen that the annual travel 
time variability along a 400 km path for the lowest 30 modes appears to be less than ±60 
ms. Additionally, for selected modes (12 through 18), the annual travel time variability 
predicted by this model was less than ±17 ms along this path. For all modes examined, 
the arrival time spectra at the 400 km range generally show a dominant frequency of 1 
cycle per year. While no attempt has been made in this investigation to separate out the 
variability due to general ocean circulation from that due to the California Current, these 
values are on the same order as that suggested by Dushaw for ray travel time variability. 
This might then suggest that the resolution achievable in this coupled ocean-acoustic 
model is not yet sufficient to accurately model the coastal region. 

The bias analysis for this environment and source placement is in agreement with 
the findings of the open ocean environment. Specifically, steeper arrivals tend to be more 
stable in the presence of mesoscale structure, as indicated by the bias results at the 400 km 
range. From bias results at the 474 km range, this environment also suggests that 
bathymetry effects may have an impact on the relative sign of the bias. Here, a relatively 
steep rise in bathymetry appeared to produce a warm shift in the bias. 
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VI. CONCLUSIONS 



The focus of this research was the bias and travel time variability of acoustic 
propagation in the presence of ocean mesoscale structure. Because such analysis is 
routinely addressed in the context of normal modes, a modal decomposition of the 
predicted acoustic fields was required. One of the major contributions of this work was 
the quantitative assessment of the errors associated with the modal decomposition of an 
approximate solution by using an improper normal mode basis set. Specifically, it was 
shown that the wide-angle PE solutions do not, in general, decompose into the standard 
normal modes of the true acoustic wave equation. Because the wide-angle PE is a 
common approximation employed in various current research efforts, and forms the basis 
of the Navy standard PE model, the development and characterization of the proper mode 
basis set is an important result for the underwater acoustics community. Although a 
completely general numerical algorithm to compute this basis set was not obtained, 
significant progress was made toward this goal. 

The main results presented in the body of this dissertation examined the bias due to 
oceanic mesoscale phenomena, and investigated its variability in terms of parameters such 
as range, local mode number, and mesoscale perturbation strength. Although previous 
investigators have attempted to address the issue of mesoscale bias using approximate 
theoretical models, those results have sometimes been ambiguous and contradictory. 
Furthermore, no model has previously been developed which provides the critical 
information on the associated variability of the bias. The deep ocean environment 
research presented here is the first to employ advanced, full-wave modelling techniques to 
explicitly compute the travel time bias introduced by mesoscale structures and its 
variability. The results suggest that nonlinear bias effects due to mesoscale ocean 
structure have a complex relationship with the sound speed perturbation field and 
bathymetry. In concluding this dissertation, a comparison of general findings from this 
work to that of earlier investigators (as discussed in Chapt. I) is provided. 
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• Previous investigations have suggested that the non-linear bias contribution is 
approximately proportional to range. Results from the open ocean 
environment with Rossby waves support this finding with the axial delay being 
on order of 50-200 ms/1000 km. 

• Flat, axial rays were found by previous investigators to be less linear than steep 
rays. In the open ocean environment analyzed here, perturbations were shown 
to produce a variable overall spreading in the arrival structure, primarily 
producing a delay in the later, axial arrivals. Earlier arrivals were shown to be 
relatively stable in the presence of mesoscale structure. 

• Earlier investigators suggested that the non-linear bias is small compared to 
seasonal fluctuations. In contrast, the analysis of a region of the California 
current system suggests that while the mean bias may be small, the bias 
variability for some environments may exceed seasonal fluctuations in travel 
time perturbations. Additionally, the degree of bias variability was found to 
have no clearly discernible relationship to range or mesoscale strength 
parameter. 

• Earlier investigations reached divergent conclusions with regard to the general 
sign of the bias. This dissertation work -suggested a predominantly cold bias 
for the open ocean environment field composed of Rossby waves. The 
acoustic path through the California Current suggested that bathymetry effects 
may have an impact on the relative sign of the bias. This then suggests a 
complex relationship between mesoscale bias and its dependence on 
bathymetry and the sound speed perturbation field. 

For the California Current model study, the primary quantities of interest were the 
modal travel times and their seasonal variability. The annual travel time variability along 
a 400 km path for the lowest 30 modes appeared to be less than ±60 ms, with some modes 
showing a variability of less than ±17 ms. For all modes examined, the arrival time 
spectra at the 400 km range generally displayed a dominant frequency of 1 cycle per year. 
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These values are on the same order as that suggested by a previous investigation using 
limited CalCOFI data for ray travel time variability. This might then suggest that the 
resolution achievable in this coupled ocean-acoustic model is not yet sufficient to 
accurately model the coastal region. 

There are several areas for potential research which logically follow this investi- 
gation. The first area deals with further development of the wide angle PE mode 
functions. Additional numerical work is required to develop proper numerical treatment 
of interfaces, bottom boundaries, and density gradients. Assessments of the mode 
functions for more complicated environments or other WAPE implementations (other than 
split-step Fourier) might also be performed. 

Finally, with respect to analysis of the region of the California Current system, the 
statistics of predicted travel time variability as determined from the coupled ocean- 
acoustic model output might be checked against measured travel time variability for 
consistency. In addition, insights from this variability study might be potentially applied 
in array configuration design for NPAL, refining tomographic techniques for this region, 
or for assessing travel time variability of ATOC transmissions through this region. 



95 



REFERENCES 



Athanassoulis, G. A., and Skarsoulis, E. K., (1995). “Arrival-time perturbations of 
broadband tomographic signals due to sound-speed disturbances. A wave-theoretic 
approach,” J. Acoust. Soc. Am. 97 (6), pp. 3575-3588. 

Chelton, D. B., and Schlax, M. G., (1996). “Global Observations of Oceanic Rossby 
Waves,” Science 272, pp. 234-238. 

Chin-Bing, S. A., King, D. B., Davis, J. A., and Evans, R. B., eds., (1993). “Lloyd’s 
Mirror - Wide Angle Propagation,” in PE Workshop II: Proceedings of the Second 
Parabolic Equation Workshop, Naval Research Laboratory NRL/BE/7181-93-0001 (US 
Government Printing Office), pp. 62-67. 

Chiu, C.-S., and Desaubies, Y., (1987). “A Planetary Wave Analysis Using the Acoustic 
and Conventional Arrays in the 1981 Ocean Tomography Experiment,” J. Phys. Oceanogr. 
17, pp. 1270-1287. 



Chiu, C.-S., Semtner, A. J., Ort, C. M., Miller, J. H., and Ehret, L. L., (1994). “A ray 
variability analysis of sound transmission from Heard Island to California,” J. Acoust. 
Soc. Am. 96 (4), pp. 2380-2388. 

Collins, M. D., (1988). “FEPE User’s Guide,” NORDA TN-365, Naval Ocean Research 
and Development Activity, Stennis Space Center, MS. 

Colosi, J. A., Flatte’, S. M., and Bracher, C., (1994). “Internal-wave effects on 1000-km 
oceanic acoustic pulse propagation: Simulation and comparison with experiment,” J. 
Acoust. Soc. Am. 96, pp. 452-468. 

DeSanto, J. A., (1977). “Relation between the solutions of the Helmholtz and parabolic 
equations for sound propagation,” J. Acoust. Soc. Am. 62 (2), pp. 295-297. 

Desaubies, Y., Chiu, C.-S., and Miller, J. H., (1986). “Acoustic mode propagation in a 
range-dependent ocean,” J. Acoust. Soc. Am. 80 (4), pp. 1148-1160. 

Dushaw, B., (1997). “The Effect of the California Current on Acoustic Transmissions 
from Pioneer Seamount,” ATOC Occasional Notes, 36, 7 January 1997. 

Feit, M. D., and Fleck, J. A., (1978). “Light propagation in graded-index fibers,” Appl. 
Opt. 17, 3990-3998. 

Garfield, N., Collins, C., Paquette, R., Rago, T., Carter, S., (1997). “Lagrangian 
Experiment of the California Undercurrent 1993-1995,” J. Phys. Oceanogr., submitted. 



97 



Godunov, S. K., and Ryabenkii, V. S., (1987). Difference Schemes, (Elsevier Science 
Publishers B. V., Amsterdam), pp. 94-109. 

Gradshteyn, I. S., and Ryzhik, I. M., (1994). Table of Integrals, Series, and Products 
(Academic Press, Inc.), pp. 966, 972. 

Hardin, R. H. and Tappert, F. D., (1973). “Applications of the split-step Fourier method to 
the numerical solution of nonlinear and variable coefficient wave equations, ” SIAM Rev. 
15, p. 423. 

Heathershaw, A. D., Stretch, C. E., and Maskel, S. J., (1991). “Coupled ocean-acoustic 
model studies of sound propagation through a front,” J. Acoust. Soc. Am. 89 (1), pp. 145- 
155. 

Holmes, E. S., and Gainey, L. A., (1991). “Software product specification for the PE 
model v.3.2,” OAML-SPS-22 (Naval Oceanography Office, Bay St. Louis, MS). 

Howe, B. M., Mercer, J. M., Spindel, R. C., Worcester, P. F., Hildebrand, J. A., Hodgkiss, 
W. S., Duda, T. F., and Flatte’, S. M., (1991). “SLICE89: A single tomography 
experiment,” in Ocean Variability and Acoustic Propagation, edited by J. Potter and A. 
Wam-Vamas (Kluwer Academic, Dordrecht), pp. 81-86. 

Jensen, F. B., Kuperman, W. A., Porter, M. B., and Schmidt, H., (1994). Computational 
Ocean Acoustics (AIP Press, New York), Chap. 5, pp. 271-341. 

Lee, D. and Botseas, G., (1982). “An implicit finite-difference (IFD) computer model for 
solving the parabolic equation,” NUSC Tech. Rep. 6659. 

Mackenzie, K. V., (1981). “Nine-term Equation for Sound Speed in the Ocean,” J. 
Acoust. Soc. Am. 70 (4), pp. 807-812. 

Mercer, J. A., and Booker, J. R., (1983). “Long-Range Propagation of Sound Through 
Oceanic Mesoscale Structures,” J. Geophys. Research, 88, Cl, pp. 689-699. 

Milder, D. M., (1969). “Ray and wave invariants for SOFAR channel propagation,” J. 
Acoust. Soc. Am. 46, pp. 1259-1263. 

Mooers, C. N. K., and Robinson, A. R., (1984). “Turbulent Jets and Eddies in the 
California Current and Inferred Cross-Shore Transports,” Science, 223, pp. 51-53. 

MODE Group, (1978). “The Mid-Ocean Dynamics Experiment,” Deep-Sea Research 25, 
pp. 859-910. 



98 



Munk, W. H., (1974). “Sound channel in an exponentially stratified ocean, with 
application to SOFAR,” J. Acoust. Soc. Am. 55, pp. 220-226. 

Munk, W. H., and Forbes, A. M. G., (1989). “Global Ocean Warming: An Acoustic 
Measure?” Am. Meteorological Soc. 19, pp. 1765-1778. 

Munk, W. H., Worcester, R, and Wunsch, C., (1995). Ocean Acoustic Tomography 
(Cambridge University Press, Cambridge), p. 88. 

Munk, W. H., and Wunsch, C., (1979). “Ocean acoustic tomography: a scheme for large 
scale monitoring,” Deep-Sea Research 26A, pp. 123-161. 

Munk, W. H., and Wunsch, C., (1983). “Ocean Acoustic Tomography: Rays and Modes,” 
Rev. of Geophys. and Space Physics, 21 (4), pp. 777-793. 

Munk, W. H., and Wunsch, C., (1985). “Biases and caustics in long-range acoustic 
tomography,” Deep-Sea Research 32 (11), pp. 1317-1346. 

Munk, W. H., and Wunsch, C., (1987). “Bias in Acoustic Travel Time Through an Ocean 
With Adiabatic Range-Dependence,” Geophys. Astrophys. Fluid Dynamics 39, pp. 1-24. 

Munk, W. H., and Zachariasen, F., (1976). “Sound propagation through a fluctuating 
stratified ocean: Theory and observation,” J. Acoust. Soc. Am. 59, pp. 818-838. 

Newhall, A. E., Lynch, J. F., Chiu, C.-S., and Daugherty, J. R., (1990). “Improvements in 
three dimensional ray tracing codes for underwater acoustics,” in Computational 
Acoustics: Ocean Acoustic Models and Supercomputing, edited by D. Lee, A. Cakmak, 
and R. Vichnevetsky (North-Holland, Amsterdam). 

Ocean Tomography Group, (1982). “A demonstration of ocean acoustic tomography,” 
Nature 299, pp. 121-125. 

Pierce, A. D., (1965). “Extension of the method of normal modes to sound propagation in 
an almost stratified medium,” J. Acoust. Soc. Am. 37, pp. 19-27. 

Porter, M. B., (1991). “The KRAKEN normal mode program,” Rep. SM-245 (SACLANT 
Undersea Research Centre, La Spezia, Italy). 

Robinson, A. R., Carman, J. C., and Glenn, S. M., (1994). “A Dynamical System for 
Acoustic Applications,” in Oceanography and Acoustics, edited by A. R. Robinson and D. 
Lee (AIP Press, NY), pp. 80-117. 

Semtner, A. J., (1995). “Modeling Ocean Circulation,” Science, 269, pp. 1379-1385. 



99 



Semtner, A.J., and Chervin, R.M. (1992). "Ocean General Circulation From A Global 
Eddy-Resolving Model," J. Geophys. Res. 97 (C4), pp. 5493-5550. 

Shang, E. C., (1989). “Ocean acoustic tomography based on adiabatic mode theory,” J. 
Acoust. Soc. Am. 85 (4), pp. 1531-1537. 

Shang, E. C., (1997). “Ocean Acoustic Thermometry and Tomography,” submitted for 
publication (University of Colorado, Boulder, CO). 

Shang, E. C., and Wang, Y. Y. (1993). “The Nonlinearity of Modal Travel Time 
Perturbation,” in Computational Acoustics - Vol. 1, edited by R. L. Lau, D. Lee and A. R. 
Robinson (Elsevier, North-Holland), pp. 385-397. 

Shang, E. C., and Wang, Y. Y., (1993). “Acoustic Travel Time Computation Based on PE 
Solution,” J. Comp. Acoust. 1 (1), pp. 91-100. 

Shang, E. C., and Wang, Y. Y., (1996). “The Impact of Mesoscale Oceanic Structure on 
Global-Scale Acoustic Propagation,” in Theoretical and Computational Acoustics edited 
by D. Lee, Y-H Pao, M.H. Schultz, and Y-C Teng, (World Scientific Publishing), pp. 409- 
432. 

Smith, A. R., Smith, K. B., and Chiu, C.-S., (1996). "Effects of Ocean Mesoscale 
Features on Adiabatic Mode Travel Times," in Proceedings of the Third European 
Conference on Underwater Acoustics, edited by J. S. Papadakis, (Crete University Press, 
Heraklion), pp. 297-302. 

Smith, K. B., (1992). “Broadband PE propagation through mesoscale perturbations,” J. 
Acoust. Soc. Am. Suppl. 91, p. 2327. 

Smith, K. B., Brown, M. G., and Tappert, F. D., (1992a). “Ray chaos in underwater 
acoustics,” J. Acoust. Soc. Am. 91, pp. 1939-1949. 

Smith, K. B., Brown, M. G., and Tappert, F. D., (1992b). “Acoustic ray chaos induced by 
mesoscale ocean structure,” J. Acoust. Soc. Am. 91, pp. 1950-1959. 

Smith, K. B., and Tappert, F. D., (1993). “UMPE: The University of Miami Parabolic 
Equation Model, Version 1.0, ” Marine Physical Laboratory Technical Memo 432. 

Spiesberger, J. L., (1985a). “Ocean acoustic tomography: Travel time biases,” J. Acoust. 
Soc. Am. 77 (1), pp. 83-100. 

Spiesberger, J. L., (1985b). “Gyre-Scale Acoustic Tomography: Biases, Iterated 
Inversions, and Numerical Methods,” J. Geophys. Research, 90 (C6), pp. 11869-11876. 



100 



Spiesberger, J. L., and Worcester, P. F., (1983). “Perturbations in travel time and ray 
geometry due to mesoscale disturbances: A comparison of exact and approximate 
calculations,” J. Acoust. Soc. Am. 74 (1), pp. 219-225. 

Spindel, R. C. and Worcester, P. F., (1990). “Ocean Acoustic Tomography,” Scientific 
American, 263 (4), pp. 94-99. 

Spindel, R. C. and Worcester, P. F., (1991). “Ocean Acoustic Tomography: A Decade of 
Development,” Sea Technology, July 1991, pp. 47-52. 

Stammer, D., Tokmakian, R., Semtner, A., and Wunsch, C., (1996). “How well does a 
1/4° global circulation model simulate large-scale oceanic observations?” J. Geophys. 
Research, 101 (CIO), pp. 25779-25811. 

Staten, R. A., Chiu, C.-S., and Semtner, A., (1996). “A Simulation Study of the 
Variability of Acoustic Transmissions from Hawaii to Monterey,” in Theoretical and 
Computational Acoustics , edited by D. Lee, Y.-H. Pao, M.H. Schultz (World Scientific 
Publishing), pp. 505-525. 

Tappert, F. D., (1977). “The parabolic approximation method, ” in Wave Propagation and 
Underwater Acoustics, Lecture Notes in Physics, edited by J. B. Keller and J. S. Papadakis 
(Springer- Verlag, Berlin), Vol. 70, pp. 224-287. 

Tappert, F. D., and Tang, X., (1995). “Acoustic travel time effects of oceanic mesoscale 
structure at global ranges,” unpublished paper, (U. of Miami, Miami, FL) dated 27 Oct 
1995. 

Thomson, D. J., (1993). “PE-Based Spectral Decomposition,” in PE Workshop II: 
Proceedings of the Second Parabolic Equation Workshop, edited by S.A. Chin-Bing, D.B. 
King, J. A. Davis, and R.B. Evans, Naval Research Laboratory NRL/BE/71 8 1-93-0001 
(US Government Printing Office), pp. 296-320. 

Thomson, D. J., and Chapman, N. R., (1983). “A wide-angle split-step algorithm for the 
parabolic equation,” J. Acoust. Soc. Am. 74, pp. 1848-1854. 

Worcester, P. and Spindel, R., (1997). “North Pacific Acoustic Laboratory,” presented at 
Office of Naval Research Long Range Acoustic Propagation Workshop, (Lake 
Arrowhead, CA), 3-4 March 1997. 



101 



INITIAL DISTRIBUTION LIST 



1 . Defense Technical Information Center 2 

8725 John J. Kingman Rd., STE 0944 

Ft. Belvoir, VA 22060-6218 

2. Dudley Knox Library 2 



Naval Postgraduate School 
411 Dyer Rd. 

Monterey, CA 93943-5101 

3. Dr. Kevin Smith, Code PH/Sk 5 

Physics Department 

833 Dyer Rd., RM 142A 
Naval Postgraduate School 
Monterey, CA 93943-5117 

4. Dr. Ching-Sang Chiu, Code OC/Ci 2 

Oceanography Department 

833 Dyer Rd.,RM 318 
Naval Postgraduate School 
Monterey, CA 93943-5122 

5. Chairman, Code PH 1 

Physics Department 

833 Dyer Rd., RM 201 
Naval Postgraduate School 
Monterey, CA 93943-5117 

6. Chairman, Code EC 1 

Electrical and Computer Engineering Department 

833 Dyer Rd., RM 437 
Naval Postgraduate School 
Monterey, C A 93943-5121 

7. Dr. Ellen Livingston, Code 3210A 1 

Office of Naval Research 

800 North Quincy Street 
Arlington, VA 22117 



103 



8. Dr. John Roese 

NCCOSC RDTE D-103 
53570 Silvergate Ave 
San Diego, CA 92152-5031 



104 



OX LIBRARY 
-TGRADUATc SCHOOL 
, ca c*mcio< 



DUDLEY KNOX LIBRARY 




